0771 32 


14690-H002-RO-00 









((M&SA-Cfi- 134272) SYSTEMS IMPSOVED 
MOHEBiCaL EIFFEBEHCISG AUaLYZEfi (SINDl) 
EUGINEEBIHG” PROGRAM MANOAJL <TRW Systems 



Group) 354 p 


00/99 




R74-75224 


Onclas 

39316 







t 


T. ISHIMOTO and L.C. FINK 




Systems ^proved Numerical Differencijig Malyzer 
Engineering-Program Manual 

T. Ishimoto and L. C. Fink 
June 1971 


NASA Contract 9-10435 



Prepared fors 

national Aeronautics and Space Administration 
Manned Spacecraft Center 
Under Contract NASA 9-10435 


Prepared by: 


T. Ishimoto 




^ 1-iJ t ^ 


L. C. Fink 


Approved by: 



R. G. Payne, Program ^Manager 


Approved by: 



R. L. Dotts 

NASA Technical Monitor 
NASA Manned Spacecraft Center 



ACKNOWLEDGEMENT 


This engineering-program manual includes the efforts of several 
individuals who in no small way contributed to the publication of this 
document. Many thanks are due to: 

J. D. Gaski , who as the prime mover behind the SINDAS 
(SINDA, CINDA-3G and CINDA, provided valuable counseling on the 
inner workings of the program, expecially the execution routines. Many of 
the computational features in these routines are original with him; it 
is through his cooperative assistance that the writing of this document 
became a reasonable venture. 

R. L. Potts, who must be given special mention not because of 
his position as NASA/MSC technical monitor but because of his sincere 
Interest in the use of SINDA and his willingness to provide assistance 
whenever and wherever possible. 

Mrs. Dorothy Gramlich, who typed this manuscript with a cheerful, 
professional attitude and technique that made this phase of the program a 
more pleasant one. 


iii 



TABLE Q? CONTENTS 


Page 

ACKNOWLEDGEMENT ill 

PREFACE xili 

IKIMENCLATURE AND MNEMONICS 1-1 

1.1 Nomenclature '* 

1.2 Mnemonics 1-2 

1.2.1 Control Constants » « 

1.2.2 ^merlcal Solution lUsutines 1-4 

1.2.3 Options 

1.2.4 Routines and Subroutines of Preprocessor 1-5 

1.2.5 Others 1-7 

2. BACKGROUND ON SIIIDA 2-1 

3. GENERAL SINDA PROGRAM DESCRIPTION 3-1 

3.1 SINDA Operating System ” 

3.2 Use of Lumped-Parameter Concept ” 

3.3 Pseudo-Compute Sequence (PCS) 3-2 

3.3.1 Long Pseudo-Compute Sequence (LPCS)^^^^^^^^^^^^^^ 3-5 

3.3.2 Short Pseudo-Gompute Sequence (SPCS) 3-6 

3.3.3 Second Long Pseudo-Compute Sequence (LPCS2) 3-8 

3.3.4 Pseudo-Compute Sequence One (PCS!) and " 

Pseudo-Compute Sequence Two (PCS2) 

3.4 Data Logistics 3-9 


3.4.1 Relative Numbers 

3.4.2 Storage Requirements and Dynamic Storage 


Allocation 

3.5 Order of Computation 3-10 

4. PREPROCESSOR 4-1 , 

4.1 General Description ” 

4.2 Description of Subroutines 4-5 

4.2.1 Routine Name: SINDA 4-6 

4.2.2 Routine Name: PREPR0 4-7 

4.2.3 Subroutine Name: ALPINT 4-8 

4.2.4 Subroutine Name; BLKCRD 4-9 

4.2.5 Subroutine Name C0DERD 4-10 


iv 






Page 

4.2.6 

Sobroutine Name; 

C0NVRT 

4-11 

4.2.7 

Subroutine Name: 

DATAED 

4-12 

4.2.8 

Subroutine Name; 

IRRHES 

4^13 

4.2.9 

Subroutine Namet 

FINDRM 

4-14 

4.2.10 

Subroutine Name: 

GENLNK 

4*15 

4.2.11 

Subroutine Name: 

GENUK 

4-16 

4.2.12 

Subroutine Name: 

INC0RE 

4-17 

4.2.13 

Subroutine Name: 

MXT0FN 

4-18 

4.2.14 

Subroutine Name: 

N0DEDA 

4-19 

4.2.15 

Subroutine Name: 

PCS2 

4-20 

4.2.16 

Subroutine Name: 

PRESUB 

4-21 

4.2.17 

Subroutine Name: 

PSEUD0 

4-22 

4.2.18 

Subroutine. Name: 

QDATA 

4-23 

4.2,19 

Subroutine Name: 

RELACT 

j 4-24 

4.2.20 

Subroutine Name: 

SEARCH 

* 4-25 

4.2.21 

Subroutine Name: 

SETFMT 

4-26 

4.2.22 

Subroutine Name: 

SINDA4 

4-27 

4,2.23 

Subroutine Name: 

SKIP 

4-28 

4.2.24 

Subroutine Name: 

SPLIT 

4-29 

4.2.25 

Subroutine Name: 

SQUEEZ 

4-30 

4,2.26 

Subroutine Name: 

STFFB 

4-31 

4.2.27 

Subroutine Name: 

TYPCHK 

4-32 

4.2.28 

Subroutine Name: 

WRTDTA 

4-33 

4.2.29 

Subroutine Name: 

WRTPMT 

4-34 

4.2.30 

Subroutine Name: 

WRTBLK 

4-35 

4.3 

Labeled Common Variables 

4-36 

4.3.1 

Labeled Common Hap 

If 

4.3.2 

Definition of Labeled Common Variables 

4-38 

4.3.3 

Dynamic Storage Structure 

4-42 

4.4 

SINDA "Tapes" and 

Their Formats 

4-46 

4.4.1 

LB3D ^ Program Data "Tape" 

It 

4.4.2 

LB4P - Program F0RTRAN "Tape" 

4-47 

4.4.3 

INTERN - Preprocessor Scratch "Tape" 

II 

4.4.4 

LUTl - Dictionary 

"Tape" 

fl 

4.4.5 

LUTS - Parameter 

Rims "Tape" 

4-48 


V 





Page 

4.5 

Overlay Structure 

4-49 

4.6 

Structure of Pseudo-Compute Sequences 

4-51 

4.6.1 

Descriptions 

ft 

4.6.2 

Structure of PCSl 

ft 

4.6.3 

Structure of PCS2 

4-53 

4.7 

Other Information 

It 

4.7.1 

Subroutine Lengths 

tt 

4.7.2 

Maximum Thermal Problem Size and Maximum 
Data Value Size 

II . . 

5. 

REVIEW OF LUMPED PARAMETER EQUATIONS AND BASIC 
NUMERICAL SOLUTIONS 

5-1 

5.1 

Lumped-Parameter Representat ion 

5-2 

5.1.1 

Some Thoughts on Lumped-Parameter Errors 

5-3 

5,1.2 

Lumped-Parameter Equations 

5-5 

5.2 

Basic Finite Difference Formulations 

5-6 

5.2.1 

Forward Finite Difference Explicit Method 

5-8 

5.2.2 

Implicit Finite Difference Method 

5-9 

5.2.3 

Steady State 

5-11 

5.2.4 

Some Comments 

fl 

6. 

SINDA NUMERICAL SOLUTION ROUTINES 

6-1 

6.1 

Objective and Presentation Format 

II 

6.2 

General Computational Procedure and Features 

6-2 

6.2.1 

Order of Computation 

tf 

6. 2. 1.1 

Finite Dif f erenee Algorithm 

ft 

6. 2. 1.2 

Updating of Optionally Specified Properties 

6-5 

6.2.2 

Operations Blocks 

If 

6.2. 2.1 

EXECUTION Operations Block 

6-9 

6.2.2. 2 

VARIABLES 1 Operations Block 

If 

«.2.2.3 

VARIABLES 2 Operations Block 

ff 

6. 2. 2. 4 

0UTPUT CALL Operat ions Block 

6-10 

6.2.3 

Control Constants 

ft 

6. 2. 3.1 

Alphabetical Listing and Brief Description of. 
Control Constants 

It 

6. 2. 3. 2 

User Specified and Optionally User . Specified 
Control Constants 

6-16 


Vi 



Page 


6.2.4 

Tlsae-Step Calculations 

6-22 

6.2.5 

Computation of Temperatures 

6-25 

6. 2.5.1 

Transient Explicit Routines 

fl 

6. 2.5.2 

Transient Implicit Routines 

6-29 

6. 2. 5.3 

Steady State Routines 

6-34 

6.2.6 

Radiation Damping 

6-37 

6.2.7 

Acceleration of Convergence by Extrapolation 
Technique 

6-39 

6. 2.7.1 

Extrapolation Technique 

II 

6 . 1 . 7. 1 

Programming Considerations 

6-41 

6 . 2 . 7. 3 

Routines Using Acceleration of Convergence 

II 

6. 2. 7. 4 

Comment on Acceleration of Convergaice 

it 

6.2.8 

Other Characteristics of the SINDA Numerical 
Solution Routines 

6-42 

6. 2. 8.1 

Units 

If 

6.2, 8.2 

General Comments on Computational Features 

II 

6.3 

Transient Explicit Solution Routines 

6-43 

6.3.1 

Subroutines: CNFRWD and CNFRDL 

6-45 

6.3. 4.1 

General Comments 

If 

6. 3. 4. 2 

Finite Difference Approximation and Computational 
Algorithm 

6-46 

6. 3.1.3 

Comments on the Computational Procedure 

6-47 

6. 3. 1.4 

Control Constants ' 

6-48 

6.3. 1.5 

Error and Other Messages 

6-49 

6.3.2 

Subroutine: CNFAST 

6-55 

6.3.2. 1 

General Comments 

If 

6. 3.2. 2 

Finite Difference Approximation and Computational 
Algorithm 

ft 

6. 3.2.3 

Comments on the Computational Procedure 

6-56 

6. 3.2.4 

Control Constants 

If 

6. 3.2. 5 

Error and Other Messages 

6-57 

6.3.3 

Subroutine: CNEXPN 

6-60 

6. 3.3.1 

General Comments 

It 

6. 3.3.2 

Finite Difference Approximation and Computational 
Algorithm 

ft 

6. 3,3. 3 

Comments on the Computational Procedure 

vii 

6-61 





Page 

6. 3. 3.4 

Control Constants 

6-62 

6. 3. 3. 5 

Error and Other Messages 

II 

6.3.4 

Subroutine; CNDUFR 

6-66 

6. 3. 4.1 

General Comments 

If 

6. 3. 4. 2 

Finite Difference Approximation and Computational 
Algorithm 

II 

6.3. 4.3 

Comments on the Computational Procedure 

6-68 

6. 3. 4.4 

Control Constants 

II 

6. 3.4. 5 

Error and Other Messages 

6-69 

6.3.5 

Subroutine; CNQUIK 

6-72 

6. 3. 5.1 

General Comments 

II 

6. 3. 5.2 

Finite Difference Approximation and Computational 
Algorithm 

II 

6. 3.5. 3 

Comments on the Computational Procedure 

6-73 

6. 3. 5. 4 

Control Constants 

6-74 

6.3. 5. 5 

Error and Other Messages 

II 

6.4 

Transient Implicit Solution Routines 

6-78 

6.4.1 

Subroutine; CNBACK 

6-80 

6. 4. 1.1 

General Comments 

II 

6.4. 1.2 

Finite Difference Approximation and Computational 
Algorithm 

II 

6. 4. 1.3 

Comments on the Computational Procedure 

II 

6.4. 1.4 

Control Constants 

6-82 

6. 4. 1.5 

Error and Other Messages 

6-83 

6.4.2 

Subrout ine : CNFWBK 

6-86 

6.4.2. 1 

General Comments 

II 

6.4. 2. 2 

Finite Difference Approximation and Computational 
Algorithm 

II 

6.4.2. 3 

Comments on the Computational Procedure 

6-88 

6.4. 2. 4 

Control Constants 

6-89 

6. 4. 2. 5 

Error and Other Messages 

II 

6.4.3 

Subroutine; CNVARB 

6-93 

6.4. 3.1 

General Comments 

II 

6. 4. 3. 2 

Finite Difference Approximation and Computational 
Algorithm 

1 

6.4. 3. 3 

Comments on the Computational Procedure 

6-95 


viii 



Page 

6. 4. 3. 4 Control Constants 6-95 

6. 4. 3. 5 Error and Other Kessages ’* 

6.5 Steady State Routines 6-98 

6.5.1 Subroutine: CINDSS 6-101 

6. 5. 1.1 General Comments •• 

6. 5. 1.2. Finite Difference Approximation and Computational 6-102 

Algorithp 

6. 5. 1.3 Comments on the Computational Procedure' 6-103 ‘ 

6.5.1.,4 Control Constants 6-104 

6.5.1.5 Error and Other Messages ” 

6.5.2 Subroutine: CIKDSL 6-107 

6. 5.2.1 General Comments ’• 

6. 5.2.2 Finite Difference Approximation and Computational " 

Algorithm 

6. 5. 2. 3 Comments on the Computational Procedure 6-109 

6. 5. 2. 4 Control Constants ” 

6. 5. 2.5 Error and Other Messages ** 

6.5.3 Subroutine: GINDSM 6-113 

6.5. 3.1 General Comments ” 

6. 5.3.2 Finite Difference Approximation and Computational " 

Algorithm 

6. 5.3. 3 Comments on the Computational Procedure 6-114 

6. 5. 3.4 Control Constants 6-115 

6. 5. 3. 5 Error and Other Messages " 

7. REFERENCES 

APPENDICES 

A.. Computer Listings of SINDA Explicit Solution 

Routines 

CNFRWD a-2 

CNFKDL A-13 

CNFAST A-24 

CNEXPN A-33 

CNDUFR A-44 

CNQUIK A-55 


lx 



Page 


Computer Listings of SINDA Implicit Solution 
Routines . 

CNBACK 

CNFWBK 

CNVARB 


Computer Listings of SINDA Steady State Solution 
Routines 

CINDSS 

CINDSL 

CINDSM 


X 



Figure No. 

3.3-1 

4.1- 1 

6.2^1 

6 . 2 - 2 

6.2- 3 

6.2- 4 ‘ ■ 

6.3- 1 

6.3- 2 

6.3- 3 

6.3- 4 

6.3- 5 

6.3- 6 

6.3- 7 

6.4- 1 

6.4- 2 

6.4- 3 

6.5- 1 

6.5- 2 

6.5- 3 


FIGURES 

Thermal Circuit for a One-Dimensional System 

SINDA Preprocessor - Major Logic and Interface 
with the User Specified Problem 

General Order of Computation for Numerical 
Solution Routines 

Computation Pattern for Numerically Solving 

Explicit and Implicit Finite Difference Algorithm 

Evaluation of Nonlinear Capacitance, Source or 
Conductance 

Method of Extrapolation to Accelerate Convergence 

QStJM and GSUM for "Block" Diffusion Node Temperature 
Calculation, CNFRWD and CNFRDL 

Calculation of Arithmetic-Node Temperature by 
"Successive Point" Iteration 

Functional Flow Chart for CNFRWD and CNFRDL 

Functional Flow Chart for CNFAST 

Functional Flow Chart for CNEXPN 

Functional Flow Chart for CNDUFR 

Functional Flow Chart for CNQUIK 

Functional Flow Chart for CNBACK 

Functional Flow Chart for CNFWBK 

Functional Flow Chart for CNVARB 

Functional Flow Chart for ClNDSS 

Functional Flow Chart for CINDSL 

Functional Flow Chart for CINDSM 


xi 


Page 

3- 4 

4- 2 

6-3 

6-4 

6-26 

6-40 

6-52 

6-53 

6-54 

6-59 

6-65 

6-71 

6-77 

6-85 

6-92 

6-98 

6-106 

6-112 

6-118 



Table Ifo, 

3.2- 1 

3.3- 1 

3.3- 2 


3.3-3 

4.1- 1 

6 . 2 - 1 
6.2-2 

6.2- 3 

6.2- 4 


6.2- 5 

6.3- 1 

6.3- 2 

6.3- 3 

6.3- 4 

6.3- 5 

6.4- 1 

6.4- 2 

6.4- 3 

6.5- 1 

6.5- 2 

6.5- 3 


TABLES 


Tabular Identification, of Conductor and Adjoining 
Node Numbers 

Example of Conductor Connections 

Long Pseudo-Compute Sequence (LPCS) for the 
Example of Table (3.3-1) 

Short Pseudo-Compute Sequence (SPCS) for the 
Example of Table (3.3-1) 

Terminology used in Description of SINDA 
Preprocessor Routines 

Optionally Specified Capacitance Expressions 

•' Impressed Source Expressions 

" ' " Coefficient Expressions 

Definition of S 3 nnbols for Tables 6.2-1 through 

6.2-3 

Characteristics of User Specified Control Constants 

Basic Computational Steps for CNFRWD and CNFRDL 

Basic Computational Steps for CNF AST 

Basic Computational Steps for CNEXPN 

Basic Computational Steps for CNDUFR 

Basic Computational Steps for CNQUIK 

Basic Computational Steps for CNBACK 

Basic Computational Steps for, CNFWBK 

Basic Computational Steps for CNVARB 

Basie Computational Steps for CINDSS 

Basic Computational Steps for CINDSL 

Basic Computational Steps for CINDSM 


Page 

3-5 

3-7 


4-2 

6-6 

tl 

6-7 

6-8 

6-11 

6-51 

6-58 

6-64 

6-70 

6-76 

6-84 

6-91 

6-97 

6-105 

6-111 

6-116 


xii 



PREFACE 


The present SINDA computer program has evolved from the CIITOA-SG 
program, which in turn evolved from the CINDA program, etc. With each 
major program revision an updated user manual was generated, but a more 
in-depth presentation of programming considerations and the theoretical 
development of the numerous subroutines were hot generated. This SINDA 
program manual represents a preliminary effort to fill some of the existing 
void by describing the pro grajn structure, by Identifying the major functions 
of each processor routine with a functional flow chart, and by a more 
in-depth mathematical description of the numerical solution subroutines. 

It is not the intent of this engineering-program manual, however, to pro- 
vide sufficient detailed information for a user to make modifications 
and/or additions to the existing subprograms. 
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NOMENCLATURE AND MNEMONICS 


1.1 


Ifomenclature 


A 

A 

(A/il) 


ij 


DD 


DN 

F,F1,F2 

"id 

k 

L 

N 

NNA 

NND 

R 

P 

% 

t 

At 

t 

m 

T 

T 

m 


= k(A/JL).. conduction coefficient between nodes i 
*3 » 

and j» 

“ array ' 

= area 

= effective ratio of cross-sectional area to distance 
between nodes i and j. 

» radiation factor between nodes i and. j (composed of 
radiation interchange factor and area) 

* capacity of ith node 

=» C^/At, capacity of ith node divided by time-step 
® 1 - DN (allows certain fraction of "old" temperature 
to be included as part of temperature change for 
current time-step, refer to Section 6. 2. 5.1) 

E DAMFA (user control constant, refer to Sections 6. 2. 5.1 
and 6 .2. 3. 2) 

= multiplying factors, either user constants or literals, 
refer to Tables 6.2-1, 6.2-2 and 6.2-3). 




(ij + T?HTj t T.) 


* a^^ or <7b^j, conduction or radiation coefficient, 

sa thermal conductivity 

** a literal multiplying factor 

= number of variable temperature nodes (NNA + NND) 

*= nianber of arithmetic-nodes 
= number of diffusion nodes 

* resistance 

= total number of nodes 
» impressed heat load into the ith node 
» time 
5* time-step 

» (TIME0 + TIMEN)/2.0, mean time 

* temperature ("F or °R) 

= (T^ + Tj)/2.0, mean temperature (“p or °R) 
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x»y 

a 

a 

e 

3’ 


» coordinate 

® (k/C), theraal diffusivity 
P 

Z G../C, , refer to equation 6.3-7 

- i 

“ factor that ranges from 0 to 1/2 (refer to equation 6.2-7) 

“ 23 (used in subroutine CNVAEB) 

» radiation interchange factor including inter- 

reflections between nodes i and j. 

—8 4 2 

“Stefan- Boltzmann constant (.1714 x 10 Btu/hr ®R ft ) 

At 
n 


At 


h-1 


, weighting factor (equation 6.3-13) 


(A^:t ) Interpolated value of array A using t as the independent variable 


m' 


m 


(aSt^) 

It 

If 

If 

If 

If 

ft 


ft 


If 

tf 

(A^;T^) 

It 

ft 

ff 

If 

It 

If 


It 



If 

(A^:T ) 

XB. 

tf 

ft 

If 

It 

ft 

ti 

T 

m 

If 



tf 


Polynomial 

ft 

ft 

If 

If 

tf 


tf 


ft 

If 


ff 

If 

If 

n 

II 

It 


It 

II 


ti 

(aP;T ) 

m 

It 

If 

If 

ft 

It 

ft 

T 

m 

It 



ft 


(A“:T^,tjj^) Interpolated value of the bivariate array A using and t^ 
as independent variables. 

(A^:T^,tjjj) Interpolated value of the bivariate array A using and t^j^ 
as independent variables . 


Subscripts 

i 

i 

ij 

i,n 

i,k 

ij.n 

ij.k 


m 


“ ith node 
= jth node 

“ between nodes i and j 

“updating of ith temperature^ source, etc. at time- 
step n. 

“ updating of ith temperature, etc. at kth iteration. 
= updating of coefficient between nodes i and j at 
time-step n 

= updating of coefficient between nodes i and J at 
kth iteration 


mean 



1.2 


l&iemouics 


1.2.1 Control constants (refer to Sections 6. 2. 3.1 and 6. 2. 3. 2) 

ARLXCA “ allowable arithmetic node relaxation temperature change 

ARLXCC = calculated maximum arithmetic node relaxation temperature change 

ATMPCA = allowable arithmetic node temperature change 

A3MPCC = calculated maximum arithmetic node temperature change 

BACKUP = back switch 

BALENG = specified system energy balance 
CSGFAC » time-step factor 
GSQdAX * maximum, value of C^/S G^^j 
CSGMIN minimum value of C./S G. . 

3. Ij 

CSGRAL « allowable range between CSQIIN and CSGMAX 
DAMPA = arithmetic node damping factor 
DAMPD *= diffusion node damping factor 

DRLXCA = allowable diffusion node relaxation temperature change 

DRLXCC = calculated diffusion node relaxation temperature change 

DTIMEH » allowable maximvim time-step 

DTIMEI * specified time-step for implicit solutions 

DTIMEL “ allowable minimum time-step 

DTIMEU = contains computed time-step 

DTMPCA 5= allowable diffusion node temperature change 

DTMPCC = calculated maximum diffusion node temperature change 

ENGBAL “ calculated system energy balance 

ITEST contains dummy integer constant 

JTEST = contains dummy integer constant 

KTEST = contains dummy integer constant 

LAXFAC = number of iterations for linearized lumped parameter system, 
CINDSM only. 

LINECT = line counter location for program output 
L00PCT = contains number of iterations performed 
N0C0PY * contains no copy switch for matrix users 
NL00P = number of specified iteration loops 
0PEITR = output each iteration switch 

PAGECT = page counter location for program output 



RTEST = contains dummy floating point constants 

StEST - contains dtirnmy floating point constants 

TIMEM » (TIKE0 + TIMEN)/2.0, mean time for computational interval 

TUdEN * TIMEN + DTIMEU, new time at the end of computational interval 

TIMEND = problem stop time 

TIME0 “ old time at the start of the computational interval 

TTEST = contains dummy floating point constants 

UTEST = contains dummy floating point constants 

VTEST = contains dummy floating point constants 

1.2.2 Mxmierical Solution Routines (refer to Sections 6.3 - 6.5) 

CINDSS = steady state routine, refer to Section 6.5.1 

CINDSL “ steady state routine, refer to Section 6.5.2 

CINDSM ® steady state routine, refer to Section 6.5.3 

CNBACK = implicit routine, refer to Section 6.4,1 

CNDUFR = explicit routine, refer to Section 6.3.4 

CNEXPN “explicit routine, refer to Section 6.3.3 
CNEAST “ explicit routine, refer to Section 6.3.2 

CNEEDL = explicit routine, refer to Section 6.3.1 

CNFWBK = implicit routine, refer to Section 6.4.2 

CNFWRD = explicit routine, refer to Section 6.3.1 

CNQUIK “ explicit routine, refer to Section 6.3.5 

CNVARB = implicit routine, refer to Section 6.4.3 

1.2.3 Options (used in Tables 6.2-1 - 6.2-3) 

BIV “ Mvariate Interpolation Variable 

BIT “ Rouble Interpolation with Time as variable 

DIV = Double Interpolation Variable 

DPV = Rouble polynomial Variable 

DTV = Double interpolation with Time and Temperature as Variables 
SIT = pingle Interpolation with Time as variable 
SIV = pingle Interpolation. Variable 
SPV “ _Single polynomial Variable 
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1.2.4 


Routines and Subroutines of Preprocessor 

SINDA = routine that specifies overlay of preprocessor to system allocator. 

PREPR0 = main routine for preprocessor; initializes counters and F0P.TRAN 
logical units; sets length of dynamic storage array and controls 
major logic. 

ALPINT = subroutine that accepts an integer in alphanumeric format and 
converts it to integer format; determines relative number of 
this actual number and converts it back to alphanumeric format. 

BLKCRD =» subroutine that formats the five generated F0RTRAN routines 

(SINDA, EXECTN, VARBU , VARBL2 , and 0UTCAL) in 507 word blocks. 

C0DERD “Subroutine that reads and checks the block header cards for the 
data blocks. 

C0NVRT « subroutine that converts Hollerith data to integer data. 

DATARD ® subroutine that scans the data block card images ijcnder an A 

\ 

format and determines appropriate format to reread the: card 
images. 

EREMES * subroutine that prints most of the error messages generated 
within the data blocks, 

FINDRM ® subroutine that moves the data in the dynamic storage array 
either up or down by 100 words ^ 

GENLNK “ subroutine that generates the driver (F0RTRAN routine named 
SINDA) for the user's program. 

GENUK » subroutine that generates user constants. 

INC0RE » subroutine that reads data into the dynamic storage array for 
the parameter-runs option. 

MXT0FN = subroutine that processes data for the "m" option (converts card 
images from mixed F0RTRAN/ SINDA notation to F0RTRAN notation. 

N0DEDA » subroutine that processes data for node and conductor data 
blocks. 

PCS2 * subroutine that packs the F0RTRAN addresses for the array and 
constants locations required by the second pseudo-compute 
sequence. 

PRESUB = subroutine that reads and checks the block header cards for the 
operations blocks and generates the non-executable F0RTRAN 
cards for each of the operations blocks via a call to BLKCRD. 



PSSUD0 » subroutine that forms the first and second pseudo-compute 
^ sequences. 

QDATA = subroutine that checks and processes all data input in the source 
data block. 

RELACT =* subroutine that finds the relative node numbers from the actual 
node niunber; computes the FORTRAN address for arrays and user 
constants from the actual number. 

SEARCH = subroutine that retains a relative number for nodes, conductors , 
user constants, and arrays, given the actual number. 

SETFMT = subroutine that processes the card for the "new format" option; 
that is, it sets up the format for data cards as specified by 
the cards with an "N" in colximn one. 

SINDA4 = subroutine that reads and processes the user input cards from 
the operations blocks. 

SKIP =» subroutine that is used when a problem is RECALLED; it positions 
the tape to the proper problem as specified on the first card of 
the data deck. 

SPLIT = subroutine that reads the data from the RECALL tape and splits 
the RECALL information onto the proper data "tape" and the 
dictionary "tape . " 

SQUEEZ = subroutine that compresses the specified data groups in the 
dynamic storage array. 

STFFB = subroutine that fills out a card image in array KBLK with 
^llerith blanks. 

TYPCHK = subroutine that Checks the input from the data blocks for the 
correct type (integer, floating point, or alphanumeric. 

WRTDTA ® subroutine that writes the program data "tape" in the format 
required by INPUTT or INPUTG. 

WRTPMT = subroutine that writes the required data for parameter runs ' 
on the parameter "runs" "tape" and the dictionary "tape." 

WRTBLK = subroutine that writes the 507 word blocks contained in array 
KBLK on the program F0RTRAN "tape." 
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Others 


SINDA = Systems Improved Numerical Differencing Analyzer 

UPCS = Long Pseudo-Compute Sequence 

SPCS * Short Pseudo -Compute Sequence 

LPCS2 =» Second Long Pseudo-Compute Sequence 

PCSl = Pseudo -Compute Sequence One 

PCS2 * Pseudo -Compute Sequence Two 

TSUM * elapsed time from -last printout 

TPRINT * time of last printout. 
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2- BACKGROUND ON SIKDA 

The original CINDA (Chrysler Improved Numerical Differencing 
Analyzer) computer program was developed by the Thermodynamics Section of 
the Aerospace Physics Branch of Chrysler Corporation Space Division at 
MSA Michoud Assembly Facility and was coded in FORTRAN-II and FAP for the 
IBM-7094 computers. CINDA was the product of an intensive analytical, 
engineering and programming effort that surveyed numerous thermal analyzer- 
type programs and studied several in-depth. The foundation for ClNDA was 
the storage and addressing of information required only for the network . 
solution and the systems features which allowed the re-utilization of core 
storage area and brought into core only those instructions necessary for 
the solution of a particular problem. A systems compiler computer program 
that automatically optimized the utilization of computer core space was 
developed. This meant the generation of an integrated operation of relative 
addressing, packing features, peripheral tape storage units and overlay 
features. 

CINDA evolved into CINDA-3G^ which was developed by the same 
group that generated CINDA with a major portion of the work done under 
contract NASA/MSC NAS9-7043. CINDA-3G represented (essentially) a com- 
plete rework of CINDA in order to take advantage of the improved systems 
software and machine speeds of the 3rd generation computers. CINDA was 
unsuitable for standard operation on third generation computers since it 
was virtually a self-contained program having its own Update, Monitor and 
Compiler. On the other hand, CINDA-3G consisted of a preprocessor 
(written in FORTRAN) which accepted the user input data and the block data 
input. The user input data was converted into advanced FORTRAN language 
subroutines and block data input was passed onto the system FORTRAN Com- 
piler. This required a double pass on data where previously only one was 

/ 

required, but the increased speed and improved software of the third 
generation machines more than compensated for the double pass. 

SINDA^ (Systems Improved Numerical Differencing Analyzer) was 
developed by the Heat Transfer and Thermod 3 mamlcs Department of TRW Systems 
Group, Redondo Beach. Most of the improvanents and subroutine additions 
to CINDA-3G was done as part of the MSA/MSC contract NAS 9-8389, 

* Superscript numbers refer to the literature cited in the Reference Section. 



entitled, "Development of Digital Computer program for Thermal Network 
Correction." Programming and systems integration were directed to the 
DNIVAC-1108 computer. 

SINDA relied quite heavily on CINDA-3G and data deck compatibility 
was rigorously followed; as a result, CINDA-3G data decks should, in the 
main, be directly operational on the SINDA program although some differences 
exist. For example, properties are updated before VARIABLES 1 call in 
CINDA-3G, whereas the properties are updated within the numerical solution 
routines after VARIABLES 1 call in SINDA. The primary differences between 
SINDA and CINDA-3G are; (1) elimination wherever possible of assembly 
language coding; (2) increased mnemonic options to aid the program user 
in data input; (3) inclusion of a second pseudo computer sequence for 
evaluation of nonlinear network elements; and (4) additional subroutines 
such as STEP (sensitivity analysis) and KAL0BS-KALFIL (Kalman filtering). 
Most of the changes and additions to CINDA-3G were required in order to 
Integrate the thermal network correction subroutine package into the I 
existing CINDA-3G program. 

During the development of SINDA a number of useful improvements 
became apparent. As a result, modifications to SINDA were made a part of 
the NASA/MSC contract NASA 9-10435 entitled, "Development of an Advanced 
SINDA Thermal Analyzer System." These changes that- include a variable 
input format, simplified parameter runs, and generated user constants 
were made by the same group that developed SINDA. These improvements are 
reported in an updated SINDA users manual."* 
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GENERAL SISDA PROGRAM DESCRIPTION 


3.1 SINDA Operatiiig System 

SINDA is more like an operating system rather than applications 
program. SINDA is programmed as a preprocessor in order to accommodate 
the desired operations relative to overlay features, data packing, 
dynamic storage allocation and subroutine library file, but yet be 
-written in F0RTRAN, This preprocessor operates in an integral fash- 
ion with a library of numerous and varied subroutines, * which may be 
called in any desired sequence but yet operate in an integrated manner. 

The preprocessor reads the input data, assigns relative numbers, packs 
this information, forms a pseudo-compute sequence(s) (which will be described 
briefly in a later paragraph of this section and is described in more 
detail in Section 4, called Preprocessor), and writes the operations 
blocks on a peripheral unit as F0RTRAN source language with all of the 
data values dimensioned exactly in labeled common. In turn, dontrols are 
shifted to the system F0RTRAN compiler which compiles the constructed 
subroutines and enters execution. The F0RTRAN allocator has access to the 
SINDA subroutine library and loads only those subroutines called by the 
problem being processed. 

As a result of this type of systems operation SINDA is extremely 
dependent upon the systems software. However, once the program is 
operational on a particular computer, the user-prepared problem data deck 
can be confined to the control cards and deck set-up requirements at a 
particular installation. 

It should be recognized that the use of a preprocessor provides 
a computer with a large capability and considerable flexibility, but 
because of the niaaerous options that are generally offered, user instruc~ 
tions are more difficult than other thermal analyzer-type programs which 
have less flexibility. 

3 . 2 Use of Lumped-Parameter Concept 

Use of SINDA is based on a lumped parameter representation of a 
physical system. This means that SINDA does not solve a set of partial 
differential equations that represents a distributive system, but rather 
SIOT)A numerically solves a set of ordinary (and in general) nonlinear 




differential equations that represent a lumped parameter system. The 
procedure for the formulation and the numerical solutions of the lumped 
parameter equations are reported extensively in literature and basic 
considerations are presented in Section 5. For the discussion to follow 
on the pseudo compute sequence it is convenient to indicate a general set 
of ordinary linear differential heat balance equations, 


dt C, 




(3.2-1) 


where. 


i = 1,2,...,N (number of variable temperatures) 
T .* constant, N < j ^ p 

= the ith nodal capacity 

q^ » the heat load into node i (impressed) . 

a^j * the conduction coefficient between nodes i and j [= ^ 


t = time 




Suppose an implicit numerical method as discussed in Section 5 .2.2 of 
this manual is chosen; the implicit finite difference form becomes after 
letting, 

dT./dt = (T. - T. )/At, T. * T. and T. * Tj 

1 x,n+l x,n ’ J 3 ,n+l x i,n+l. 


At ^'^i,n+l ■* ^i,n^ ” ^i ^ij ,n+l " ^i,n+l^ 

where, n “ temperature at time point t^ 

^i n+1 ~ tCTiperature at time point t^^_^ 

At = time-step 

Rearrangement of equation (3.2-2) yields, 

_ P N _ p 

(C. + I a,.) T, . - - Z a.. T. = q. + C. T. + Z 
i ij i»n+l xj j,n+l ^i x x,n 

j?^i j?*i 

i » 1,2,...,N 


(3.2-2) 


a. , T, 
ia j ,n 


(3.2-3) 


T. = constant, N < j < p 
J 

where, = C^/At, average capacity of node 1 over At time-step 

3*3 Pseudo-Compute Sequence (PCS) 

A pseudo-compute sequence as generated by the SINDA. preprocessor 



is a list of numbers that indicates the position of required data values 
in various arrays such as conductance, temperature and capacitance. This 
meaning will become clearer by formulating equation (3.2-3) in a matrix 
form. The matrix formulation is straightforward since temperatures at 
time-step n+1 are the unknowns and terms on the right side of equation 
(3.2-3) represent the forcing function. Let us expand equation (3.2-3) to 
show this. 


_ P 


(C,+ E a, .)T. a,., T- j^,..,-a,„ T„ * q,+ G, T_ + E a,. T. 

1 Ij 12 2, n+1’ IN N,n+1 U 1 l,n J 


P 

E 

j=N+l 


_ P _ P 

"*a*- T- , . + (C/^+ E a«, . )T— , - > * * T., , • — qj.+ C* T— + E a* , T , 

21 l,n+l 2 2j' 2, n+1* 2N N,n+1 ^2 2 2,n 2j j ,n 


"®N1 ^1 ,n+l"®N2 ^2 ,n+l » * * ’ ^Sl’*' . ^ ®Nj ^"^NjU+l 

3“x 

Thus the matrix form of equation (3.2-3) becomes, 

m {T»} = {Q} 


C., T-, + E a... T. 

^ J— N+1^^^ 3 


where, 


g = 


j«2 


-a, 


21 


-a, 


N1 


‘IN 


^2j^*'’*’ ’■®2N 

in • : 

P _ 

» “®N2 ****V-f^% 

5n 


(3.3-1) 


(3.3-2) 


l,n+l 


q, + G, T. + E a-. T. 


- ' ’^ 2 .n+l \ ». 3 - 3 ) J {Q} -j ^ + Cj I ( 3 . 3 - 4 ) 


N,n+1 


q.., + G„ T„ + E a-.. T. 

■ N N,n “Nj J,n 


The matrix represented by equation (3.3-2) appears to, be a full 
matrix (very small number of elements that are zero), but in reality most 
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of the off-diagoaal elements are zero. Thus, if equation (3.2-3) was to 
be solved by a matrix Inversion technique, all elements including zeros 
must be stored. Since the number of elements varies as (N is the number 
of nodes), the required number of data locations would Vary as and the 
computer time required for matrix inversion would be proportional to N^. 

The explicit and iterative implicit numerical methods (refer to 
Section 5) of solving equation (3.2-1) lend themselves for optimizing the 
data storage area required and for reducing the solution tine. If the 
conductors are numbered and related to the appropriate adjoining nodes as 
Indicated in Table (3.2—1), retention of adjoining node number for each 
conductor provides a means of identifying element position in the coefficient 
matrix. This can be seen by considering the one-dimensional heat conduc- 
tion example pictured in Figure (3.3-1). 


I 


IW m fit fp •1’ 

1 Gx2 2 G23 5 



X ^ X 


Figure 3. 3-1-. Thermal Circuit for a One-Dimensional System 


The set of equations -associated with the problem of Figure (3.3-1) 
may be readily expressed as. 


(Cj+Gi). _ -Gi , 0 , 0 , 

t (, 02 ^ 2 ^ 2 ^ * ~^2 * ® » 

0 , ^2 -Gy , 

0 , 0 , -63 ,(C^-W34G^), 


(3.3-5) 


» "^4 » *^5^4) 


By cou^arlng the element position of equation (3.3-5) with the 
tabular identification in Table (3.2-1) it is seen that elements with zero 
values need not be stored. The main diagonal term is never zero and is a 
composite of capacitance and off-diagonal conductors. 
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Table (3.2-1) Tabular Identification of Conductor 
and Adjoining Node Numbers 


Conductance 

ith 

Adjoining 


Number 

Node 

Node Number 

Comment 

G# 

N# 

N# 



G1 is qonductor #1 between 
nodes 1 and 2. 

G1 is conductor #1 between 
nodes 2 and 1. . 

G2 is conductor #2 between 
nodes 2 and 3. 

G4 is conductor #4 between 
nodes 5 and 4. 

It is of interest to note that the use of a pseudo-compute sequence 

■! 

is only one of a number of ways to store data efficiently ^ For example, 

TRW TAP® does not employ a pseudo-compute sequence because of other user 
requirements. However, from a data storage standpoint, it appears that 
the use of a pseudo-compute sequence utilizes computer core most efficiently. 

More than one pseudo-compute sequence is formed by SINPA. Both a 
so-called long (LPCS) and a so-called short (SPCS) pseudo-compute sequence 
as used in CINDA-3G^ are formed and in addition a second long pseudo-compute 
(LPCS2) required for thermal network correction is also formed in SINDA. A 
detailed discussion of these pseudo-compute sequences will be presented in Sec- 
tion 4.6, but is of interest here to indicate the characteristics of these 
"sequences^” 

3.3.1 Long Pseudo-Compute Sequence (LPCS) 

A long pseudo-compute sequence identifies the position and 
value of all off-diagonal elements of the coefficient matrix. This is done 
by operating on adjoining node nijanbers which have been assigned relative 
node numbers by the preprocessor. Since nodal temperatures are calculated 
sequentially in ascending numerical order, the conductor and adjoining node number 
are searched until node one is found with the conductor ntunber and the 
other adjoining node number stored in a single core location. In addition, 
several indicators are stored in this single core location. These 
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1 1 2 

1 2 1 

2 2 3 

• • • 

• • • 

4 5 4 


< 



i Indicators are : (1) yar C (Indicates the input of a capacitor as a vari- 

able); (2) var G (indicates the input of a conductor as a variable); (3) rad 
(indicates the input of a radiation conductanee) ; (4) Q (indicates the 
input of. a source in the source data block); (5) one-way (indicates the input 
of a one-way conductor); and (6) last G (indicates the last conductor to a 
particular node). Order of indicator storage is indicated in Table (3.3-2). 

Search is continued imtil all node-one's have been located and 
characteristics processed. The procedure is repeated for all node-two's 
and sc forth sequentially until all nodes have been processed. The 
important consideration of a LPCS is the encounter of each conductor of 
the coefficient matrix twice. Formation of a pseudo-compute sequence for 
the sample shown in Table (3*3-1) is given in Table (3.3-2). A pseudo- 
compute sequence starts with node one and advances the node number by one 
each time a last conductor Indicator (last G) is passed. The conductor 
and node numbers identify the position of the conductor value in an array 
of conductor values and the position of the temperature, capacitor and 
source values in arrays of temperature, capacitor and source values 
respectively. 

A long pseudo-compute sequence is well-suited for "successive 
point" iteration (refer to Section 5.2.2 for a discussion of this) of the 
implicit finite difference, equations because all elements of the coefficient 
matrix are identified. Thus, when a row of the coefficient matrix is 
processed and a new value of temperature obtained, the new temperature 
can then be used in the calculation procedure of succeeding rows. 

3.3.2 Short Pseudo-Compute Sequence (SPCS) 

The short pseudo-compute sequence identifies each conductor only 
once and since the coefficient matrix (equation 3.3-1) is symmetrical, all 
sparsity and off-diagonal elements of the coefficient matrix are accounted 
for. The node being processed and the adjoining node nximber reveal 
temperature- and source-value locations. The short pseudo-compute sequence 
for the example in Table (3.3-1) is formed in Table (3.3-3). By placing 
a minus sign on the initially encountered other-adjoining nodes, these nodes 
are not recognized on a second encoimter. A short pseudo-compute sequence 



Table (3.3-1) Example of Conductor Connections 


Conductor No . 
G# 

1 

2 

3 

4 

5 

6 


Adjoining Node Numbers 


N# 

1 

1 

1 

2 

2 

3 


N# 

2 

3 

4 

3 

4 
4 


Table (3.3-2) Long Pseudo-Compute Sequence (LPCS) 
for the Example of Table (3.3-1) 


Node No. 

Last 

var 

var 

rad Q 

G# 

One- 

Node # 

Searched 

G 

C 

G 



way 

Stored 

1 





1 


2 

1 





2 


3 

1 

1 




3 


4 

2 





1 


1 

2 





4 


3 

2 

1 




5 


4 

3 





2 


1 

3 





4 


2 

3 

1 




6 


4 

4 





5 


2 

4 

1 




6 


3 


Table 

(3.3-3) 

Short 

Pseudo-Compute Sequence 

(SPCS) 




for 

the Example of Table 

(3.3-1) 



Node No. 

Last 

var 

var 

rad Q 

G# 

One- 

Node # 

Searched 

G 

G 

G 



way 

Stored 

1 





1 

• 

2 

1 





2 


3 

1 

1 




3 


4 

2 





4 


3 

2 

1 




5 


4 

3 

1 




6 


4 

4 

1 




0 


0 


3 - 7 






is well-suited for explicit numerical solutions methods which calculate the 
energy flow through the conductor, add it to the source location of the node 
being processed and subtract it from the source location for the adjoining 
node. The SPCS can be used for implicit methods of solution but the "block" 
Iterative procedure (refer to Section 5.2.2 for a discussion of this) must be 
used since succeeding rows of conductor and adjoining node numbers do not 
contain the necessary element information, 

3.3.3 Second Loug Pseudo-Compute Sequence (LPCS2) 

The second long pseudo-compute sequence (LPGS2) as a user input 
option flags a non-linear conductor between two diffusion nodes twice; 

LFCS flags the non-linear conductor only one. LPCS2 is required for the 
thermal network correction of a sparse network by the use of subroutine 
KAFIL (refer to Reference 3 or 6). 

3.3.4 Pseudo-Compute Sequence One (PCS 1) and Pseudo-Compute Sequence 

Two (PCS 2) 

PCS 1 and PCS 2 are not user options but are fixed internally. 

The contents of PCS 1 and PCS 2 are governed by the user input of LPCS, 

SPCS or LPCS2). PCS 1 contains two relative addresses (conductor and 
adjoining node locations), two non-linear type indicators, and an Impressed 
source indicator. Indicators are keyed through a simple counter to a 
second pseudo-compute sequence (PCS 2) which contains integer addresses 
or relative constant and array starting locations necessary for evaluation 
of temperature varying coefficients and time varying coefficients for 
sources. When the input data contain literal values in SIV type calls, 
the preprocessor stores the values as extended user constants and 
supplies the relative constant address to the second pseudo-compute 
sequence. Detailed discussion on PCS 1 and PCS 2 is presented in Section 4.6. 
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Data Logistics 
3«4f.l Relative Numbers 

Both the long and short pseudo-compute sequences require the 
storage of only the finite values in the coefficient matrix, thereby 
talcing advantage of matrix sparsity. If the short pseudo-compute sequence 
is used, the advantage of sjTnmetry is accounted for. Gonductors with the 
same constant value may share the same conductor number and value. The 
storage efficiency of the pseudo-compute sequences requires the sequential 
numbering of the nodes and the conductors. Since the numbering of thermal 
math-models is arbitrary and not sequential, the SINDA program assigns 
relative numbers . (starting from one, sequential and ascending) to the 
actual numbers of the incoming node data, conductor data , constants data 
and array data in the order received . Thus, numbers not used in the 
actual numbering system are neither identified nor required. 

3.4,2 Storage Requirements and Dynamic Storage Allocation 

All numerical solution subroutines require three locations for 
each diffusion node data (tenperature, capacitance and source), two 
locations for each arithmetic node data (temperature and source), one 
location for each boundary data (temperature) and one location for each 
conductor value, In addition intermediate data storage ranging from zero 
to three locations per node may be required for the storage of temperatures 
and temperature differences; acceleration of convergence (refer to 
Section 6.2.7) used in the implicit and steady state routines (except 
GINDSS) requires three locations. Storage requirements for conductances 
depends upon the problem. For example, each internal diffusion and 
arithmetic node of a three-dimensional conduction system with rectangular- 
nodalization will be connected with only three being unique; thus, each 
diffusion node (or arithmetic node) in a three-dimensional conduction system 
requires from six to nine storage locations for data values (temperature, 
capacitance, source j three conductors and up to three intermediate locations). 
Now each of the conductors for the short pseudo-compute sequence requires 
a single core location that contains two integer values (conductor and 
adjoining node numbers) and six Indicators (refer to Section 3.3.1 for 
description). Each of the conductors between variable temperatures for the 



long pse«do-compute sequence requires two core locations since the con- 
ductors are used twice during the computational process. This means that 
each internal node of a three-dimensional conduction system will require 
six data addressing locations for the long pseudo-compute sequence and, 
on the average, three data addressing locations for the short pseudo- 
compute sequence. 

Thus for a three-dimensional conduction system (no radiation), 
the nximber of required core locations per node can vary ffom nine (tempera- 
ture, capacitance, source, three unique conductors and three data addressing 
locations) to fifteen (temperature, capacitance, source, six conductors 
and six data addressing locations) exclusive of the second pseudo-compute 
sequence which is required for variable coefficients, capacitance and 
sources. 

The user must allocate an array of data locations which is to be 
used for intermediate data storage and initialize the array start and length 

indicators. Each subroutine that requires intermediate storage area has 

* 

access to this array and the start and length indicators. During a sub- 
routine execution a check on the sufficiency of space is made and start 
and length indicator are updated. If a subroutine calls upon another 
subroutine that requires intermediate storage, the called Subroutine 
repeats the check and update procedure. Whenever any subroutine terminates 
its operation, the start and length indicators are returned to their entry 
values. This process is termed "Dynamic Storage Allocation" and allows 
subroutines to share a common working area. 

3.5 Order of Computation 

A network data deck consists of four data blocks (node, conductor, 
constants, and array), one optional data block (source) and four operations 
blocks which are preprocessed by the preprocessor and passed on to the system 
F0RTRAN compiler. Non-network problems require no node or conductor data 
blocks. The operations blocks are named EXECUTI0N, VARIABLES 1, VARIABLES 2 
and OUTPUT CALLS; the SINDA preprocessor constructs these blocks into 
individvial subroutines with the entry names EXECTN, VAEBLl, VARBL2 and 
0UTCAL, respectively. After a successful F0RTRAN compilation, control is 
passed to the EXECTN subroutine. This means that the order of computation 
depends on the sequence of subroutine calls placed in the EXECUTI0N block 
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by the program user. No other operations blocks are performed unless 
called upon by the user either directly by name or indirectly through a 
subroutine call. The numerical solution subroutines described in 
Section 6 internally call upon VARBLl, VARBL2 and 0UTCAL; The Internal 
order of computation for these routines is similar with the primary 
difference being the numerical solution method. A general flow diagram 
of the numerical solution routines, as well as a detailed description of 
each is presented in Section 6. 
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4. 


PREPROCESSOR 


4.1 General Description 

The SINDA preprocessor reads and fflcialyzes the user input deck 
and from this information constructs a program tailored to the user*s 
requirements. 

The rationale for a preprocessor is flexibility and speed. 
Flexibility is achieved by providing the user with a library of routines 
to solve problems, manipulate data, and print selected values. In addition, 
the user may insert non-SIKDA routines into the constructed program. Speed 
(defined here as minimal execution time) is achieved by structuring the 
data in an efficient manner. 

The SINDA preprocessor consists of thirty routines with seven 
overlay links. All of the routines are written in F0RTRAN except for one 
assembly language routine which writes a "tape" in a format acceptable to 

I 

the F0RTRAN compiler. These routines provide the user with a nxjmber of 
major options in the type of problem to be solved and the form of the 
data to be used. Henceforth these major options are designated as "major 
logic" of the preprocessor. See Figure 4.1-1 for a flow chart of the 
major logic of the preprocessor and its interface with the user program. 

The major logic consists of the five following options: (1) NASA 

MSG EDIT feature; (2) RECALL option; (3) generation of a THERMAL 
problem; (4) generation of a GENERAL problem; (5) and PARAMETER RUNS 
option. The primary features of each item of the major logic is discussed 
below. 

(1) EDIT feature: The first card of the deck is checked for 

the user request of the EDIT feature. If the EDIT feature 
is requested the input "tape" is changed from the system 
input "tape" to the EDIT "tape" and control is transferred 
to subroutine EDIT for processing. On return a branch is 
made to the THERMAL or GENERAL section as specified by the 
data on the EDIT "tape. " If the EDIT feature has not been 

requested, the check for RECALL is made. 

RECALL option: The first' card of the deck is checked for 

user request of the RECALL option. If the RECALL option 
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Read & edit data; 
write "tape" for 
input to SINDA 


Initialize dynamic storage 
length and F0RT1AN log ical 
units 


Recall 


Find proper problem 
read data from "tape" 


Read data blocks, write data 
& dictionaries on "tape" 


GENERAL 


THERMAL 


Form pseudo-compute sequence, 
write it on "tape" 


Compile the five 
F0RTRAN routines 
from "tape" & 
any load & go 
subroutines 


Write users main program 
"SIlffiA" on "tape” 


Read instructions blocks, write 
four corresponding subroutines 
on tape 


Read data and 
dictionaries 
from "tape" 


Execute the 

user’s 

program 


Process parameter runs, if any, 
write data on "tape" 


End of 
S^ata ^ 


Terminate preprocessor 


Figure 4.1-1. SINDA Preprocessor - Major Logic 
and Interface with the User Specified Problem 
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is requested, control is transferred to subroutine SPLIT for pro- 
cessing. On return a branch is made to the PAEAMETER RUNS section. 
If the RECALL option has not been requested, the second card of 
the deck is checked for the type of problem, THERMAL or GENERAL. 

(3) THERMAL problem; The type of pseudo-compute sequence requested is 
noted, the title block is read, the data blocks are read and pro- 
cessed, the pseudo-compute sequence is formed, the driver for the 
user program (SINDA) is written on "tape," the operations blocks 
are read and processed and their F0RTRAN equivalents are written 

■on "tape," and finally a check is made for the user requests of 
the PARAMETER RUNS option. 

(4) GENERAL problem: This section is identical to a THERMAL problem 

except that only constants data and array data of the data blocks 
are read and processed; a pseudo-eompute sequence is not 
formed. 

(5) PARAMETER RUNS option: A check is made for the user request of 

the PARAMETER RUIJS option. If the PARAMETER RUNS option is 
requested, the appropriate data blocks are read and processed. ■ 

If not, the preprocessor is terminated. 

Description of SINDA preprocessor routines is presented in the 
sections to follow. Terminology used in the description is listed and 
defined in Table 4.1-1. 



Table (4.1-1) Terminology Used in Description of 
SINDA Preprocessor Routines 

(1) DATA BLOCKS: The five user input blocks which contain data rather 

than instructions; these DATA BLOCKS are NODE DATA, CONDUCTOR DATA, 
CONSTANTS DATA, ARRAY DATA and the optional blocks SOURCE DATA. 

(2) OPERATIONS BLOCKS; The four user input blocks which contain instruc-' 
tions on problem solution, as opposed to data contained in the DATA 
BLOCKS. These OPERATIONS BLOCKS are EXECUTI0N, VARIABLES 1, 

VARIABLES 2 and 0UTPUT CALLS. 

(3) Non-fatal error: An error that does not terminate the preprocessor 

Immediately. That is, the preprocessor will continue scanning the 
remaining cards of the input deck for errors. However, the user 
program will not be executed. 

(4) Fatal error; An error that terminates the run Immediately. 

(5) N/A; Means not ^plicable. 

(6) "TAPE"; The term "tape" in quotes is used to signify any external 
storage device. That is, any piece of computer hardware, excluding 
the central processor, on which data can be stored and retrieved. The 
three most familiar examples are; magnetic tape, drian and disk. 

(7) Dictionary; A list of the actual SINDA manbers in relative order. 

For example, the actual node number corresponding to the kth 
relative node number is the kth item of the node nianber dictionary. 

(8) Data group: A data group composed of the pertinent Information 

extracted from a particular data block. For example, the two groups 
derived from the constants data are; the user constants numbers and 
the user constants values. 

(9) Bit manipulation; Terminology that implies the ability to store and 
access information within a computer word. This capability is also 
called packing and vinpacking. 

(10) Routine; A general term used to describe any program element. 

(11) Subroutine: A special type of program element that is callable from 

a routine. 

(12) Fixed constants; The term used in the preprocessor for control constants. 
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4.2 


Description of Subrotjttines 

Sections 4.2.1 through 4.2.30 below describe the 30 routines 
of the SINDA preprocessor. The descriptions are based on the UNIVAC 1108 
computer under the EXEG II operating system; it should be understood, 
however, that much of the information is machine-dependent and is depen- 
dent upon the facilities operating system. Note that the element named 
SINDA (Section 4.2.1) and the element n^ed PREPR0 (Section 4.2,2) are not 
subroutines in the technical sense of the word; hence, these two elements 
are referred to by the more general term ’’routine." 

Each element of the preprocessor is described by the following* ‘ 
eleven subtitles; 

(1) ” SUBROUTINE NAME - this specifies the name of the element. 

(2) PROGRAMMING LANGUAGE - This may be F0RTRAN, ASM, or MAP. F0RTRAN 

implies F0RTRAN V, ASM stands for assembly language (some- 
times called SLEUTH II) and MAP is a special language which 
defines the overlay structure, 

(3) PURPOSE - This gives a brief statement of the factional 

capabilities of the element. 

(4) RESTRICTIONS - This gives an indication of where the input 

parameters come from, the form of the input parameters and 
the placement of the output parameters. 

(5) "TAPES" USED - This represents a list of each F0RTRAN logical 

unit referenced within this element.- 

(6) SPECIAL FEATURES — This specifies programming features that are 

unique to a particular machine. 

(7) OTHER SUBROUTINES CALLED -This represents a list of the 

external references. : 

(8) CALLING SEQUENCE - This gives a list of the subroutine argviments, 

if any, and a brief discussion of their use. 

(9) ERROR PROCEDURES - This discusses the steps taken when an error 

is encountered. 

(10) STORAGE REQUIRED - This gives the octal and decimal storage 

required for this element. 

(11) LABELED COMMON - This represents a list of each labeled common 

name used in this element. 
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4.2.1 


ROUTINE NAME: SINDA 


PROGRAMMING LANGUAGE ; MAP 

PURPOSE ; Tliis routine specifies the overlay structure of the preprocessor 
to the system allocator (loader). 

RESTRICTIONS ; N/A 

”TAPES" USED ; N/A 

SPECIAL FEATURES ; N/A ^ ' 

OTHER SUBROUTINES USED ; N/A 
CALLING SEQUENCE ; N/A 
ERROR PROCEDURES ; N/A 
STORAGE REQUIRED ; N/A 


LABELED COMMON; N/A 



4.2.2 


ROUTINE NAME; PEEPR0 


PROGRAMHING LANGUAGE ; F0RTRM 

PURPOSE ; This routine is the main routine (i.e., the driver) for the pre- 
processor. It initializes the counters and FORTBAN logical units, sets the 
length of the djmamic storage array, and controls the major logic. The 
major logic includes; (1) the EDIT feature (NASA MSG only); (2) the 
RECALL of a stored problem; (3) setup of a new user problem; (4) and 
preprocessor termination procedures. 

RESTRICTIONS ; N/A 

"TAPES'* USED; 


System input "tape" 

NIN 

System output "tape" 

N0UT 

Problem data "tape" 

LB3D 

Problem FORTRAN "tape" 

LB4P 

Dictionary "tape" 

LUTl 

Parameter runs "tape" 

LUT3 

Recall "tape" 

LUT7 

Internal scratch "tape" 

INTERN 


SPECIAL FEATURES ; System error termination - the problem data unit (LB3D) 
and the problem FORTRAN unit (LB4P) are flagged to stop before the data 
scan begins in the event that a system error terminates the preprocessor 
prematurely. The reason the problem date unit is flagged to stop is that 
for a RECALL problem the problem FORTRAN unit must not be written on. 

OTHER SUBROUTINES USED : C0DERD, GENLNK, PRESUB, PSEUD0, SINDA4, SPLIT and 

HRTBLK. 

CALLING SEQUENCE ; N/A 

ERROR PROCEDURES; The error termination procedures are controlled by three 
flags named ERDATA, PR0GRM, and ENDRUN. The three flags are ^ the labeled 
common block named DATA. ERDATA is used to flag non-fatal errors encountered 
while reading the data blocks, while PR0GEM performs the same function for 
the operations blocks. See Section 4.7.2* 

STORAGE REQUIRED ; 443 octal words = 291 decimal words. See Section 4.7.1. 
lABELED COMMON; BUCKET, CRDBLK, DATA, L0GIG, PL0GIC, and TAPE, 
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4.2.3 


SUBROUT33gE NAME; ALPINT 


PDRPOSE ; This subroutine accepts an integer in the alphanumeric format nAl, 
and converts it to integer format, determines the relative number of this 
acttial number, and converts the relative number back to an alphanumeric 
format of the form mAl, 


RESTRICTIONS ; The input and output is transmitted via the labeled common 
block named CIMAGE(see Section 4.3). The input must consist exclusively 
of the ten decimal digits. - 


*'TAPES" USED ; System output "tape” 

SPECIAL FEATURES ; None 

OTHER SUBROUTINES USED; SEARCH 


N0UT 


CALLING SEQUENCE ; 

KLET 


ALPINT (KLET, I ST , lEND , J) 

is an integer variable that indicates which dictionary 
is to be used for converting actual to relative, 
is the starting location of the alphanumeric integer, 
is the ending location of the alphanumeric integer, 
points to the last location + 1 of the converted Integer. 

In the event that a given actual number has no relative 
number in the dictionary list, an error message will be issued and the 
relative operations blocks error flag (PR0GRM) will be set to 1.0. 


1ST 

lEND 

J 

ERROR PROCEDURES; 


STORAGE REQUIRED ; 665 octal words * 437 decimal words. See Section 4.7.1. 
LABELED COMMON ; BUCKET, CIMAGE, DATA P0INT, and TAPE. 
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4.2.4 


SUBROUTINE NAME: BLKCBD 


PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE : This subroutine formats the five generated F0RTEAN routines 

(SINDA, EXECTN, VARBLl, VAR3L2, and 0UTCAL) in 507 word blocks, which are 
acceptable to the F0RTRAN compiler . This information is stored in labeled 
common block CRDBLK, array KBLK. A complete discussion of the required tape 
fo'rmat is found in UNIVAC 1108, EXEG II, Programmers Reference Manual, 
UP-4058 C, Appendix D.4 entitled. Program Elements on Magnetic Tape (via 
CUR). 

RESTRICTIONS ; The input is Hollerith card images with a 14A6 format. 

It is transmitted either through the array IMAGE in labeled common CRDBLK, 
or through "tape" INTERN. 

"TAPES" USED ; Internal scratch "tape" INTERN 
SPECIAJL FEATURES ; None 

OTHER SUBROUTINES USED: STFFB and WRTBLK 

CALLING SEQUENCE ; BLKCBD 
ERROR PROCEDURES ; none 

STORAGE REQUIRED ; 753 octal Words *491 decimal words. See Section 4.7.1. 
LABELED COMMON; GBDBLK and TAPE, 
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4.2.5 


SUBROUTINE NAME; C0DEED 


PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE : This subroutine reads and checks the block header cards for the 

data blocks. It also performs the following functions: (1) the second 

data card of the deck is checked for a thermal or general problem, and if 
it is a thermal problem the type of pseudo-compute sequence specified is 
noted; (2) the title block is read and processed; (3) the actual array 
and constant numbers from the automated options are converted into F0RTRAN 
addresses; and (4) the parameter run block header cards are read and 
checked. 

RESTRICTIONS : None 

'*TAPES USED : System input "tape" NIN 

System output '-tape" N0UT- 

F0RTRAN V reread 30 j 

SPECIAL FEATURES : The F0RTRAN V intrinsic function FLD is used for bit 

manipulation. 

GALLING SEQUENCE : C0DERD 

ERROR PROCEDURES ; In general, the errors checked for in this subroutine 
are of the fatal type; for example^ data blocks out of order. The result 
of a fatal error is that the fatal error flag (ENDRUN) is set to 1.0 and 
control is returned to PREPRO for immediate termination. 

STORAGE REQUIRED ; 3213 octal words * 1675 words decimal. See Section 

4.7.1. 

LABELED COMMON ; BUCKET, DATA, L0GIC,, PL0GIC, and TAPE. 
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4.2.6 SUBROPTINE NAME ; G01iVRT 
PROGRAMMING LAHGUAGE ; F0RTRAN ^ ' 

PURPOSE ; This subroutine converts Hollerith data to integer data. 


RESTRICTIONS ; The Hollerith data taust be contained in one word and con- 
sist of only the ten decimal digits 

”TAPES” USED ; None 

SPECIAL FEATURES ; The F0RTRAN V intrinsic function FLD is used for bit 
manipulation. 


OTHER SUBROUTINES USED; EREMES 


CALLING SEQUENCE ; C0NVRT(IST,IEND,ITEMP,CRDERR) 

1ST is the pointer to the first bit of the first character. 

lEND is the pointer to the first bit of the last character. 

IT01P is the word containing the Hollerith data on entry and 
the Integer number on return. 

CRDERR is a logical error flag which is set true if an error 
is encountered during the conversion. 


ERROR PROCEDURES ; If a non-integer is encountered, an error message is 
printed and CEDERR is set to true. 


STORAGE REQUIRED ; 150 octal words = 104 decimal words. See Section 4.7.1. 
LABELED C014M0N; None 
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4.2.7 


SUBROUTINE NAME; DATARD 


PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine scans the data block card images under an A 
format and determines the appropriate format (of the form Fn, In, or An) 
to reread the card image. The card Images are then reread under the 
generated format. In addition, the constants data block and the array 
data block are processed, 

RESTRICTIONS ; None 

“TAPES” USED ; System input "tape” NIN 

System output "tape” N0UT 

F0RTRAN V reread 30 

SPECIAL FEATURES ; The FORTRAN V intrinsic function FLD is used for bit 
manipulation. j 

OTHER SUBROUTIIifES USED ; ERRMES , FINDRM, GENUK, N0DEDA (and its entry 
point C0NDDA), SETFMT, SQUEEZ, and TYPCHK. 

CALLING SEQUENCE ; UATARD 

ERROR PROCEDURES : All errors checked for in this subroutine are non-fatal. 

An error message is printed either internally or from subroutine ERRMES and 
the data blocks error flag (ERDATA) is set to 1,0, 

STORAGE REQUIRED ; 5344 octal words ** 2788 decimal words. See Section 4.7.1. 

LABELED COMMON : BUCKET, CHECKD, DATA, FLAGS, L0GIG, PL0GIC, P0INT, and 

TAPE. 
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4.2.8 


SUBROUTINE NAME: ERBMES 


PROGRAMMING LANGUAGE; F0RTRAN 


PURPOSE ; This subroutine prints most of the error messages that can be 
generated within the data blocks. 


RESTRICTIONS ; None 

“TAPES" USED ; System output "tape" 

SPECIAL FEATURES; None ’ 


N0UT 


OTHER SUBROUTINES USED ; System subroutine EXIT 
CALLING SEQUENCE ; ERRMES (JUMP, I,J,K) 

JUMP is an integer that points to the appropriate error message 


I 

J 

K 


via a computed G0 T0 statement 

are data words which allow a maximum of three printed 
data words per error message. t 


ERROR PROCEDURES ; If the number of error messages printed exceeds 199, 
the preprocessor is terminated by a call to EXIT. 

STORAGE REQUIRED ; 2166 octal words * 1142 decimal words. See Section 4.7.1. 
LABELED COMMON; DATA and TAPE ^ ^ ' 
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4.2.9 SUBROUTINE NAME : FINDRM 

PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine moves the data in the dynamic storage array 
either up or down by 100 words. In the process it may delete certain 
groups of data that are no longer needed. 

RESTRICTIONS ; None 

"TAPES" USED ; System output "tape" N0UT 
SPECIAL FEATURES ; None 

OTHER SUBROUTINES USED ; SQUEEZ, and system subroutine EXIT. 

CALLING SEQUENCE FINDRM (L0CN0,M) 

L0CN0 is a pointer to a portion of the dynamic storage array 
where the data group that needs more room resides. 

M is the address where the next data value is to be 

stored. 


ERROR PROCEDURES ; If the dynamic storage array is full, an error message 
is printed and the preprocessor is terminated via CALL EXIT. 

STORAGE REQUIRED ! 407 octal words * 263 decimal words. See Section 4.7.1. 

LABELED COMMON; BUCKET L0GIC P0INT, and TAPE. 


4 - 14 


34 < 



4.2,10 SUBROUTI!^ NAME ; GEHLNK 
PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine generates the driver, FORTRAN routine name and 
SINDA, for the user's program. 

RESTRICTIONS ; None 

"TAPES" USED ; Internal Scratch"tape" INTERN 

SPECIAL FEATURES ; None 

OTHER SUBROUTINES USED ; BLKCED 

CALLING SEQUENCE ; GENLNK 

ERROR PROCEDURES ; None 

LABELED COmON ; CRDBLK, BATA, L0GIC, PL0Gie, and TAPE. 


4 - 15 

3S< 



4,2.11 


SUBROUTIHE HAME: GENUK 


PROGRAMMING LANGUAGE : F0RTRAN 

PDRPOSE ; This subroutine is used to goierate user constants. 

RESTRICTIONS : The input data is taken from array TiMP in labeled common 

CHECKD and the output data (i.e. , the gaoLerated user constants) is put 
into array B in labeled common BUCKET. 

**TAPES” USED ; None 

SPECIAL FEATURES ; None ' 

OTHER SUBROUTINES USED ; ERSMES, FINDEM, and TYPCHK. 

CALLING SEQUENCE ; GENUK(IWRDS) 

IWRDS is the number of words to be processed in array TEMP. 

ERROR PROCEDURES ; The input data is checked and if an error is foimd, 

control is transferred to subroutine ERRMES. 

STORAGE REQUIRED ; 451 octal words =297 decimal words. See Section 4.7.1. 
LABELED COMMON ; BUCKET, CHECKD, /DATA, and P0INT. 
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4.2.12 SUBROUTINE HAHE ; iaC0RE 
PROGRAmiNG LANGUAGE ; F0RTEAN 

PURPOSE ; This subroutine reads data into the dynamic storage array for 
the parameter rims option. 

RESTRICTIONS: None 

"TAPES" USED ; Dictionary "tape" - LUTl 
Parameter runs "tape" LUT3 

SPECIAL FEATURES ; None 

OTHER SUBROUTINES USED ; None ' 

CALLING SEQUENCE ; INC0RE (ITEST) 

ITEST is an integer flag which determines the data group 
group to be read. 

ERROR PROCEDURES ; None 

STORAGE REQUIRED ; 600 octal Words ==384 decimal words. See Section 4.7,1 
LABELED COMMON ; BUCKET, DATA, L0GIC,PL0GIC, P0INT, and TAPE. 
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) 4.2.13 SUBROUTIHE NAME ; MXT0FN 
PROGRAMMING LMGUAGE ; F0RTRAN 

PURPOSE ; This subroutine processes the data for the "M" option. That is, 
it converts card images from mixed F0RTRM/SINDA notation to F0RTEAN 
notation. 

RESTRICTIONS ; The input (array IM0LL) and output (array JH0LL) are both 
in labeled common CIMAGE and they are both in an 80A1 format. The F0RTRAN 
from array JH0LL is copied to array IMAGE tinder a 14A6 format for pro- 
cessing by the F0ETR^ compiler. 

“TAPES” used ; None 

SPECIAL FEATURES ; The FORTRAN V intrinsic function FLD is used for bit 
manipulation. , 

OTHER SUBROUTINES USED ; ALPINT and BLKCRD 

CALLING SEQUENCE ; MXT0FN 

I ERROR PROCEDURES ; None 

STORAGE REQUIRED ; 522 octal words — 338 decimal words. See Section 4.7.1. 
LABELED COMMON; CIMAGE and CRDBLK 
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4.2.14 


SUBROUTINE NAME; N0DEDA 


PROGRAMING LANGUAGE F0RTEAN 

PURPOSE ; This subroutlae processes the data for the node and conductor 
data blocks. 

RESTRICTIONS ; The input is received via labeled common CHECKD: array TEMP 
and the processed data are stored In the dynamic storage array. 

“TAPES” USED ; None 

SPECIAL FEATURES ; This subroutine has a second entry point named C0NDDA. 
Also, the F0RTRAN V intrinsic function FLD is used for bit manipulation. 

OTHER SUBROUTINES USED ; ERRMES, FINDRM, RELAGT, and TYPCHK 

GALLING SEQUENCE ; N0DEDA(JUMP,lWRDS) 
or C0NDDA(JUMP,IWRDS) 


JUMP is a flag which indicates which code option (columns 8, 
9, and 10 of the data card) the user has selected. 

IWRDS is the ntanber of data values in array TEMP to be 
processed. 

ERROR PROCEDURES ; If an error is detected while scanning the input data, 
control is transferred to subroutine ERRMES, . 

STORAGE REQUIRED ; 7030 octal words 3608 decimal words. See Section 4.7.1. 
LABELED COMtiON ; BUCKET, CHECKD, DATA, FLAGS, and P0INT. 
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4.2.15 


SUBROUTINE NAME: PCS2 


PROGRAMMING LANGUAGE : F0RTRAN 

PURPOSE ; This suhroutine packs the F0RTRAN addresses for the array and 
constants locations required by the pseudo-compute sequence. 

P.SSTRICTIONS ; None 

♦'TAPES*' USED ; None 

SPECIAL FEATURES ; The F0RTRAN V intrinsic function FLD is used for bit 
manipulation. 

OTHER SUBROUTINES' USED; None 

PCS2(IB,IPCS,LITA) 

is the word in the dynamic storage array where the 
addresses are found. 

is the word into which the addresses are packed, 
is a flag that is set to 1 if the array address was 
input as a literal and therefore has been added to 
the constants data. 

None 

54 octal words = 44 decimal words. See Section 4.7.1. 
None 


CALLING SEQUENCE ; 
IB 

IPGS 

LITA 

ERROR PROCEDURES ; 
STORAGE REQUIRED ; 
LABELED COMMON; 
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4.2.16 SUBROUTINE NAME ; PRESUB 
PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine reads and checks the block header cards for the 
operations blocks and generates the non-executable F0RTRAN cards for each 
of the operations blocks via a call to BLKCRD. 

RESTRICTIONS ; None 

•'TAPES” USED ; System input, "tape" NIN 

System output "tape" N0UT 

SPECIAL' FEATURES ; None 

OTHER SUBROUTmES USED ; BLKCRD^^^^ 

CALLING SEQUENCE ; PRESUB(N) 

N is an integer from 1 to 4 which Indicates which opera- 

tions block is being processed, 

ERROR PROCEDURES ; If the card read is not the correct block header card, 
an error message is printed and the fatal error flag is set to 1.0. 

STORAGE REQUIRED ; 200 octal words ^ 128 decimal words. See Section 4.7.1 

LABELED COMMON ; CRDBLK, DATA, L0GIC, and TAPE. 



4.2.17 SUBROUTINE NAME; PSEUD0 


PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine forms the first and second pseudo-compute 
sequence. See Section 4.6. 

RESTRICTIONS ; The necessary input is extracted from the dynamic storage 
array and the output (the two pseudo-compute sequences) is placed in the 
dynamic storage array. 

”TAPES" USED ; N0RT 

SPECIAL FEATURES ; The F0RTRAN V intrinsic function FLD is used for bit 
manipulation. 

OTHER SUBROUTINES USED ; FINDEM, PCS2 and WRTDTA. 

CALLING SEQUENCE ; PSEUD0 

ERROR PROCEDURES ; If an error is encountered while forming the pseudo- | 
compute sequences} an error message will be printed and the non— fatal 
error flag (ERDATA) is set to 1.0. 

STORAGE REQUIRED ; 2244 octal words =1188 decimal words. See Section 4.7. 

BUCKET, DATA, L0GIC, PL0GIC, and TAPE. 


LABELED COMMON; 



4.2.1S SPBROUTINE MME; QDATA 
PROGRAMMING LMGUAGE ; F0RIEAN 

PURPOSE ; This subroutine checks and processes all data input in the source 
data block. 

RESTRICTIONS ; The input is received from the calling sequence and labeled 
common CHECKD, The processed data is placed in the dynamic storage array. 

”TAPES” USED ; None 

SPECIAL FEATURES ; The F0RTRAN V intrinsic function FLD is used for bit 
manipulation. 

OTHER SUBROUTINES USED ; ERRMES, FINDRM, RELACT, and TYPCHK. 

CALLING SEQUENCE ; QDATA(C0DE,IMDS) 

C0DE is the three letter option from columns 8, 9, and 10 
of the data card. | 

IWRDS is the number of words in array TEMP to be processed. 

ERROR PROCEDURES ; If an error is encountered, control is transferred to 
subroutine ERRMES . 

STORAGE REQUIRED ; 1655 octal words = 941 decimal words. See Section 4.7.1 
LABELED COMMON ; BUCKET, CHECKD, and P0INT. 



4.2.19 


SUBROUTINE NAI-IE: RELACT 


PROGRAMMING LANGUAGE ; F0RTEAN 

PURPOSE ; This subroutine finds the relative node nxmiber from the actual 
node nzmber. In addition, it computes the F0RTEAN address for arrays and 
user constants from the actual number. 

RESTRICTIONS ; This subroutine is used in conjunction with the data 
blocks. 

"TAPES" USED ; None 

SPECIAL FEATURES ; The F0RTEAN V intrinsic function FLU is used for bit 
manipulation. 

OTHER SUBROUTINE USED ; ERRlfflS 

CALLING SEQUENCE ; RELACT (K, MM, J,JJ) 

K determines the path through the program via a computed 
G0 T0 statement. 

MM is the actual number on entry, and the F0RTRAN address 
on return. 

J ) are print variables for subroutine ERRMES. 

JJ ^ 

ERROR PROCEDURES ; In the event an error is encountered, control is 
transferred to subroutine ERRMES. 

STORAGE REQUIRED ; 311 octal words 201 decimal words. See Section 4.7.1. 
LABELED COMMON ; BUCKET, DATA, and P0INT. 
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SUBROUTINE NAME; SEARCH 


PROGRAMMING LANGUAGE; F0RTRAN 

PURPOSE; This subroutine returns a relative number for nodes, conductors, 
user constants, and arrays, given the actual nxmber. 

RESTRICTIONS ; This subroutine is used in conjxmction with the operations 
blocks. 

"TAPES" USED ; None 

SPECIAL FEATURES ; The F0RTRAN V Intrinsic function FLD is used for bit 
manipulation. 

OTHER SUmOUTINES USED ; None 

CALLING SEQUENCE ; SEARCH (N,IA,NDIM,L0C) 

N is the actual number. 

lA is the first word of the dictionary of actual numbers 
to be searched. 

NDIM is the number of words of lA to be searched. 

L0C is the relative ninnber returned to the calling program. 

ERROR PROCEDURES ; If the input actual ntimber is not found in the 
dictionary, L0C is set to zero. 

STORAGE REQUIRED ; 101 octal words = 65 decimal words. See Section 4.7.1. 
LABELED COMMON; None 
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4.2.21 SUBROUTINE NAME ; SETEMT 
\ • . ' 

' PROGRAMMING LANGUAGE ; F0RTRAN ^ 

PURPOSE ; This subroutine processes the cards for the "new format" option; 
that is» it sets up the format for data cards as specified by the cards 
with a N in colimm one. 

RESTRICTIONS ; The input/output array is passed through the calling 
sequence argument. 

"TAPES" USEP ; F0RTRAN V reread 30 

SPECIAL FEATURES ; None 

OTHER SUBROUTINES USED ; None : 

CALLING SEQUENCE ; SETFl® ( JW , B) 

JUMP is an integer flag that determines the path through 

the code. j 

B is an array which contains the card images. 

ERROR PROCEDURES ; None 

^ STORAGE REQUIRED ; 221 octal words - 145 decimal words. See Section 4.7.1. 
LABELED COMMON; None 
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4.2.22 


SUBROUTINE NAME: SINDA4 


PROGRAMMING LANGUAGE : F0RTRAN 

PURPOSE : This subroutine reads and processes the user fiiput cards froia the 

operations blocks. 

RESTRICTIONS : None 

"TAPES" USED : System input "tape" MIN 

System output - "tape" N0UT 

Internal scratch "tape" INTERN 
F0RTRAN V reread 30 

SPECIAL FEATURES : None 

OTHER SUBROUTINES USED : BLKCRD, MXT0FN, and SEARCH. 

CALLING SEQUENCE : SINDA4 (NAME) 

NAME is an integer flag that tells the subroutine! which 
operations block is being processed. 

ERROR PROCEDURES : In the event an error is encountered while processing 

the operations blocks, an error message is printed and the error flag 
PROGRM is set to 1.0. 

STORAGE REQUIRED : 2372 octal words = 1274 decimal words. See Section 4.7.1. 

LABELED COMMON : BUCKET, CIMAGE, CRDBLK, DATA, L0GIC, PL0GIC, P0INT, and 

TAPE. 
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4.2.23 


SUBROIJTIKE NAME ; SKIP 
PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine is used when a problem is being RECALLed, It 
positions the tape to the proper problem as specified on the first card 
of the data deck. 

RESTRICTIONS t The data is read from tape R. There is no output. 

”TAPES” USED ; RECAIX "tape" LUT7 
SPECIAL FEATURES ; None 
OTHER SUBROUTINES USED ; None 
CALLING SEQUENCE ; SKIP 
ERROR PROCEDURES ; None 

STORAGE REQUIRED : 324 octal words “212 decimal words. See Section 4.7.1. 

LABELED COMMON; TAPE ^ 



4.2.24 


SUBROUTINE NAME; SPLIT 


PROGRAMMING LANGUAGE ; - 

PURPOSE: This subroutine reads the data from the RECALL tape and splits 

the RECALL information onto the program data "tape'* (LB3U) smd the 
dictionary "tape" (LUTl) . 


RESTRICTIONS : The input is from the RECALL "tape" and the output is 

placed on the program data *'tape," the diet lonairy "tape," and the 
parameter runs "tape." ' 

"TAPES" USED ; RECALL "tape" ^ LUT7 

Program data "tape" LB3D 

Dictionary "tape^' LUTl 

Parameter rros "tape" . LUT3 

SPECIAL FEATURES ; None 

OTHER SUBROUTINES CALLED ; SKIP 

CALLING SEQUENCE ; SPLIT(ID) 

ID is the RECALL name punched in the first card of the 

data deck. 


ERROR PROCEDURES ; None 

STORAGE REQUIRED ; 746 octal words - 486 decimal words. See Section 4.7.1. 
LABELED COMMON ; BUCKET, DATA, and TAPE. 
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4.2.25 SUBROUTINE NAME ; SQUEEZ 
) EBOGRAMMING LANGUAGE ; F0RTEAN 

PURPOSE ; This subroutine compresses the specified data groups in the 
dynamic storage array. The compression is accomplished by placing the 
data groups sequentially in the dynamic storage array. 

RESTRICTIONS; None 


"TAPES" USED: None 

SPECIAL FEATURES: 

None 

OTHER SUBROUTINES 

USED: None 

CALLING SEQUENCE; 

SQUEEZ (1ST, lEND) 

1ST 

is the data group number where the compression is to 
start . 

lEND 

is the last data group number for this compression. 

ERROR PROCEDURES: 

None 

STORAGE REQUIRED; 

115 octal words = 77 decimal words. See Section 4.7.1. 

LABELED COMMON: 

BUCKET and P0INT 
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4.2.26 


SUBROUTINE NAME: STEFB 


PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine fills out a card image in arr^ KBLK with 
Hollerith blanks. 

RESTRICTIONS ; The pointers to the words to set to blank are in the 
calling sequence, and the array containing the card images is in labeled 
conmon CRDBLK. 

“TAPES" USED ; None 

SPECIAE FEATURES ; None 

OTHER SUBROUTINES USED ; None 

CALLING SEQUENCE ; STFFB(I,J) 

I is the first work in KBLK to set to blank. 

J is the last word in KBLK to set to blank. 

EPJROR PROCEDURES ; None 

STORAGE REQUIRED ; 41 octal words =! 33 decimal words. See Section 4.7.1. 
LABELED COMMON; CRDBLK. 
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A. 2, 27 


SUBROUTINE NAME ; TYPCHK 
PROGRAM^^ING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine checks the input from the data blocks for the 
correct type; type means integer, floating point, or alphanumeric. Also, 
it regulates the conversion of the A's and K’s for the automated options 
via a call to C0NVRT. 

RESTRICTIONS ; The input and output are transferred through the calling 
sequence arguments and labeled common CHECKD. 

•’TAPES” USED ; None 

SPECIAL FEATURES ; The F0RTRAN V intrinsic ftmctlon FID is used for bit 
manipulation. 

OTHER SUBROUTINES USED ; C0NVRT and ERRMES 

CALLING SEQUENCE ; TYPCHK(JUMP,IERR, J) 

JUHP indicates what type the word should be. 

XERR tells subroutine ERRMES which error message to print 

if the word is not of the type indicated by JUMP. 

J is a pointer to the word type in array KFLFX. 

ERROR PROCEDURES ; If a word is not of the proper type control is 

transferred to subroutine ERRMES to print an error message and the 
logical flag CRBERR is set to true. 

STORAGE REQUIRED; 233 octal words « 155 decimal words. See Section 4.7.1. 
LABELED COMMON; CHECKD 
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4.2.28 SUBROUTINE KMIE ; NRtDTA 
PROGRAMMING lANGUAGE : F0RTRAN 

PURPOSE ; This subroutine writes the program data "tape" in the format 
required by INPUTT or INPUTG, 

RESTRICTIONS ; The data to be written on "tape" is found in the dynamic 
storage array, 

"TAPES” USED ; Program datd "tape" LB3D 

SPECIAL FEATURES ; None 

OTHER SUBROUTINES USED ; None 

CALLING SEQUENCE ; mTDTA(JUMP)^^ — 

JUMP is an integer flag that indicates which data block to 
write and what format to use. 

ERROR PROCEDURES ; None 

STORAGE REQUIRED ; 645 octal words 421 decimal words. See Section 4.7.1. 
LABELED COMMON ; BUCKET, DATA, L0GIC, PL0GIC, P0INT, and TAPE 
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4.2.29 


SUBROUTINE NAME ; NRTPMT 
PROGRAMMING LANGUAGE ; F0RTRAN 

PURPOSE ; This subroutine writes the data that is needed for parameter 
runs on the parameter rims "tape" and writes the dictionary "tape." 

RESTRICTIONS ; The information that is written on the "tapes" is found in 
the dynamic storage array, 

"TAPES" USED ; Dictionary "tape" . LUTl 
Parameter runs "tape" LUT3 

SPECIAL FEATURES ; None 

OTHER SUBROUTINES USED ; None 

CALLING SEQUENCE ; WRTPMT(JUMP) ; 

JUMP is an integer flag that indicates to NRTPMT which set 
of information to write. 

ERROR PROCEDURES ; None 

STORAGE REQUIRED ; 401 octal words = 257 decimal words. See Section 4.7.1 
LABELED COMMON; BUCKET, DATA, L0GIC, P0INT, and TAPE. 
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4.2.30 SUBROUTINE NAME ; V3RTBLK 

PROGRAMMING LANGUAGE : Assembly Language 

PURPOSE ; This subroutine writes the 507 word blocks 
KBLK on the program F0RTRAN "tape." 

RESTRICTIONS ; None 

"TAPES" USED ; Program F0RTRAN "tapei" LB4P 

SPECIAL FEATURES ; None , ^ 

OTHER SUBROUTINES USED ; None 
CALLING SEQUENCE ; WRTBLK 
ERROR PROCEDURES ; None 

STORAGE REQUIRED ; 14 octal words =12 decimal words 
LABELED COMMON; CRDBLK 


I 
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contained in array 


See Section 4.7.1. 



4.3 


LABELED C0MM0N VARIABLES 


The SIRDA preprocessor uses nine labeled coimnon blocks to pass 
data and flags between the various subroutines. Labeled common names, 
in alphabetical order, are: 

BUCKET CHECKED CIMAGE 

CRDBLK DATA FLAGS 

L0GIG PL0GIC P0INT 

Note that the UNIVAC 1108 version does not utilize blank common. 
The two sections that follow give: 1) a map of the labeled common usage by 

subroutine name and by overlay link; 2) a definition of the variables 
used within each labeled common block; and 3) dynamic storage structure. 

4.3.1 Labeled Common Map 

The map below gives the labeled common name, a list of the over- 
lay links that use it by link number and a list of the routines that use it. 


LABELED C0MM0N 

0VERLAY 

R0UTINE 

NAME 

LINK NAMES 


NAMES 

BUCKET 

0, 1, 2, 4, 5 

ALPINT 

C0DERD 



DATARD 

FINDRM 



GENUK 

INC0RE 



N0DEDA 

PREPR0 



PSEUD0 

QDATA 



RELACT 

SINDA4 



SPLIT 

SQUEEZ 



WRTDTA 

WRTPMT 

CHECKD 

1 

DATARD 

GENUK 



N0DEDA 

QDATA 



TYPCHK 


CIMAGE 

4 

ALPINT 

MXT0FN 



SINDA4 


CRDBLK 

0, 3, 4 

BLKCRD 

GENLNK 



MXT0FN 

PREPR0 



PRESUB 

SINDA4 



STFFB 

WRTBLK 
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LABELED C0MM0N 

0VERLAY 

R0UTINE 

NAME 

LINK NAMES 

NAMES 

DATA 

0, 1, 2, 3, 4, 5 

ALPINT 

C0DERD 



DATARD 

ERRMES 



GENLINK 

GENUK 



INC0RE 

N0DEDA 



PREPR0 

PRESUB 



PSEUD0 

RELACT 



SINDA4 

SPLIT 

- 


WRTDTA 

WRTPMT 

fiLags 

1 

DATARD 

N0DEDA 

L0GIC 

0. 1. 2 

C0DERD 

DATARD 



FINDEM 

INC0RE 



PMPR0 

PSEUD0 



WRTDTA 

WRTPMT 

PL0GIC 

0, 1. 3, 4 

C0DERD 

DATARD 



GENLNK 

INC0RE 



PREPR0 

SINDA4 



WRTDTA 


P0INT 

0, 1, 2. 4 

ALPINT 

C0DERD 



DATARD 

FINDRM 



GENUK 

INC0RE 



N0DEDA 

PREPR0 



PSEUD0 

QDATA 



RELACT 

SINDA4 



SQUEEZ 

WRTDTA 



WRTPMT 


TAPE 

0. 1, 2. 3, 4, 5 

ALPINT 

BLKCRD 



C0DERD 

DATARD 



ERRMES 

FINDRM 



GENLNK 

INC0RE 



PREPR0 

PRESUB 



PSEUD0 

SINDA4 



SKIP 

SPLIT 



WRTDTA 

WRTPMT 
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Definition of Labeled Common Variables 

(1) Labeled common name BUCKET. 

BUCKET is the dynamic storage array (see Section 4.3.3). 

(2) Labeled common name CHECKD 

CEECKD is used to temporarily store and check the user’s 
input data. 

VARIABLE NAME DESCRIPTION 

TEMP (35) A. temporary storage array, which contains 

or ITEMP the user’s input data as read from the 

data cards. 

XGEN(35) A temporary storage array used to store a 

or IGEN copy of TEMP when the user is generating 

data. 

KFLFX(35) An indicator used to check the data array 

which contains one of the following 
numbers: 

1 * floating point number 
0 ® integer number 
'-1 » Hollerith word 

CRDERR A logical flag: Set true if and only if 

an error was found on the data card now 
being processed. 

(3) Labeled common name CIMAGE 

CIMAGE is used to store and manipulate Hollerith card images 
for the ’M' option. 

VARIABLE NAME DESCRIPTION 

IH0LL(8O) The input card image, read under an 80A1 

format . 

JH0LL(I6O) The constructed card image, also an 80A1 

format. 
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(4) Labeled common name CBDBLK 

CEDBLK is used to construct the five generated F0RTRA2J 
routines. 


VARIABLE NAME 


DESCRIPTION 


LSTART 


A logical flag that signals the start 
of n new routine if set to true. 


LECARD 


LC0PY 


A logical flag that signals the end of 
a routine. 

A logical flag that tells the program 
to copy the card image (14A6) found in 
IMAGE to the next available slot in KBLK. 


NW 


KBLK (507) 


Is a counter whose value is the next 
available word in KBLK. 

An array that contains F0RTRAN card 
images of the generated routines. 


IMAGE (14) An array that contains one card image to 

be copied into KBLK. 

(5) Labeled common name DATA 

DATA is used to store the counters that Indicate (to the 
program) how many of each data type has been encotintered. 
In addition, it contains three error flags. 


VARIABLE NAME 

NND 

NNA 

NNB 

NNT 

NGL 

MGR 

NGT 

NtJC 

NECl 


DESCRIPTION 

The number of diffusion nodes. 

The number of arithmetic nodes. 

The ntmaber of boundary nodes. 

The total number of nodes. 

The nuinber of linear conductors. 

The number of radiation conductors. 

The total number of conductors. 

The number of user constants. 

The number of added constants from 
automated options in the node data 
block. 
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NEC2 

NCT 

LENA 

ERDATA 

PR06RAM 

ENDRUN 

LSEQl 

LSEQ2 

L0NG 


The number of added constants from the 
automated options in the conductor data 
block. 

The total number of constants. 

The total number of words used in the array 
data block. 

The non-fatal error flag ,for the data 
blocks. ERDATA 0.0 means an error 
has been found. 

The non-fatal error flag for the opera- 
tions blocks. FR06RAN 0.0 indicates an 
error condition. 

The fatal error flag for the preprocessor. 
ENDRUN 5 * 0.0 signals the program to terminate 
immediately. 

The length of the first pseudo-compute , 
sequence. I 

The length of the second pseudo-compute 
sequence. 

A logical flag set to true if the user 
is requesting the long pseudo-compute 
sequence. 


(6) Labeled common name FLAGS 

FLAGS contains three flags that are used to go to the proper 

block of coding in subroutine N0DEDA. 

VARIABLE NAME ‘ DESCRIPTION 

LEAP Used with the GEN option. 

N0NLIN Flags a set of multiply connected con- 

ductors as radiation} if set to true. 

INDX Determines path when multiply connected 

conductors require more than one data 
card. 


(7) Labeled common name L0GIG 

L0GIC contains a number of logical flags and the fifty 
fixed constants. 
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VARIABLE NAME 

DESCRIPTION 

LN0DE 

Set to true if any node data was processed . 

LC0ND 

Set to true if any conductor data 
was processed. 

LC0NST 

Set to true if any user constants were 
processed 

LAERAY 

Set to true if any array data was processed. 

LPRINT 

Debug print flag, set to true if there is 
an asterisk in coltimn 80 of the BCD 3 . 
THEBMAL/GENEEAL card. 

KBRNCH 

An integer that specifies which data 
block is being processed. 

Fixe (50) 
IFIXC 

The array that contains the fixed 
(control) constants. 

KTPKNT 

Optional print flag for a list of relative 
versus actual user constant numbers. Set 
time if there is an asterisk in column 80 
of the BCD 3C0NSTANTS DATA card. 

AYPENT 

Optional print flag for a list of actual 
array numbers versus F0RTRAN address. Set 
true if there is an asterisk in column 80 - 
of the BCD 3ARRAY DATA card. 

GENEBL 

Set true for a general problem. 

LQ 

Set true if any data was processed from the 
source data block. 


(8) Labeled common name FL06IC 

PL0GIC contains a number of logical flags that are used in 
conjtinction with parameter runs. 


VARIABLE NAME 

PARINT 

PABFIN 

PN0DE 

PC0I® 

PC0NST 


DESCRIPTION 

Set true for initial parameters run. 

Set true for final parameters run. 

Set true if node data was processed. 

Set true if conductor data was processed. 

Set true if user constants data was 
processed. 
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PAERAY 


Set true if array data was processed. 

• PTITLE Set true If a new title was input, 

PCHGID Contains the alphanumeric word INITIAL 

or FINAL to be used as the rtm identification 
on ”tape^' LB3D. 

(9) Labeled common name P0INT 

P0INT is used in conjunction with dynamic storage array 
BUCKET. See Section 4.3,3, 

4.3.3 Dynamic Storage Structure 

Dynamic storage represents one of the techniques of maximizing 
problem size with a computer with finite core. In dynamic storage each 
data, set is placed sequentially into one array end-to-end. This 
eliminates the wasted core inherent with the traditional system of 
dimensioning each variable at some fixed length. However, the price 
paid for the additional core is the extra time required to compute the 
address of a variable. 

The SINDA preprocessor used three arrays to store and address 
the data sets. The data sets are stored in an array named B, or IB, or BB. 
This array resides in labeled common BUCKET. The length at which B can be 
dimensioned depends; on the system that the computer facility uses. At 
NASA MSC approximatEly 30,000 words are allocated to B, In addition, in 
labeled common P0INT there are two arrays named L0C and LEN, each 
dimensioned at 20. L0C (I) contains the starting location in B for the 
Ith data set and LEN (I) contains the length of the Ith data set. 

The infomnation below gives, in detail, the contents of the 
dynamic storage array for each data block as it exists just after the data 
block has been processed. 

(1) Node data block 
data set 1: 

bit 1, automated option flag 

bit 2, Q from S0URCE DATA flag 

bits 16-35, actual node number 
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xlata set 2: 


bits 0-35, 

temperature value 

data set 3: 


bits 0-35, 

capacitance value 

data set 4: 


bits 0-5 

non-linear capacitance t 3 n?e 

bit 6, 

literal array flag 

bits 7-20, 

actxial array number 

bit 21, 

literal constant flag 

bits 22-35, 

actual constant number 

data set 5: 


bits 0-35, 

literals encountered in 4. 

Source data block 


data set 2 (first 

word of group): 

bits 0-5, 

source option type 

bits 6-20, 

relative node number 

bits 21-35, 

not used 

data set 2 (second word of group) : 

bits 0-5, 

not used 

bit 6, 

literal array flag 

bits 7-20, 

actual array ntnnber 

bit 21, 

literal constant flag 

bits 22-35, 

actual constant number 

data set 3: 


bits 0-35, 

literals encountered in 2 

Conductor data block 

data set 6: 


bits 0-35, 

actual conductor number ' 

data set 7: 


. bit 0, 

amilti connections flag 

bit 1, 

radiation flag 

bit 2, 

automated option flag 
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(4) 


(5) 


( 6 ) 


bits 3-5, 

not used 

bit 6, 

1 way flag for NA 

bits 7, 20, 

relative node nuniber NA 

bit 21 

1 way flag for NB 

bits 22-35, 

relative node number NB 

data set 8: 

bits 0-35, 

conductance value 

data set 9: 


bits 0-5, 

conductor option type 

' bit 6, 

literal array flag 

bits 7-20, 

actual array number 

bit 21, 

literal constant flag 

bits 22-35, 

actual constant number 

data set 10: 

bits 0-35, 

literals encountered in 9 

Constants data block 


data set 11: 


bits 0—35, 

actual constant number 

data set 12: 


bits 0-35, 

constant value 

Arr^ data block 

data set 13: 


bits 0-35, 

actual array number 

data set 14: 

bits 0-35, 

array length 

data set 15: 


bits 0-35, 

array value 

Pseudo-compute sequences 

data set 16 (1st pseudo-compute sequence): 

bit 0, 

last conductor flag 

bit 1, 

automated capacitance flag 

bit 2, 

automated conductance flag 


4 - 44 



bit 3, 
bit 4, 
bits 5-20, 
bit 21 
bits 22-35, 


radiation conductance flag 
Q from source block flag 
relative conductor ntxmber 
1 way conductor flag 
relative adjoining node number 


data set 17 (second pseudo-compute sequence); 
bits 0-4, automated option type 

bit 5, not used 

bits 6-21, F0RTRAN address for array 

bit 22, not used 

bits 23-35, relative constant number 

The bit numbering convention above conforms to the UNIVAC 
standard notation, where each 36 bit word is numbered 0 through 35 from 
left to right. Each of the 1 bit flags above is querried in the following 
manner: 0 means NO, and 1 means TES. If the literal array flag or the 

literal constant flag is set to 1, then the bits immediately to the right 
of the flag do not contain the actual array or constant number. Instead, 
they contain a pointer to the next data set where the literal value is 
stored. In those data sets that store information for the automated 
options it is sometimes necessary to use more than one word per option. 
When this is the case, the automated option type (bits 0-5) is set to 0. 
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4.4 


SINDA "Tapes** and Their Formats 

The SINDA program in its normal operating mode utilizes six 
"tapes." Five of these "tapes" are assigned by the program and the re- 
maining one contains the program; it is assigned via control cards. The 
store and recall options require one additional "tape" each and the NASA 
edit feature requires two additional "tapes." The following paragraphs 
contain information on the five normal SINDA "tapes," 

^ *4 . 1 LB3D - Program Data "Tape" ‘ 

This "tape" is set up by the preprocessor (WRTDTA) and read 
by INPUTT, for a thermal problem, or INPUTG, for a general problem, just 
prior to performing the instiructions of the execution block. The contents 
of this unit are: 

(1) Problem identification. 

WRITE (LB3D)RUNID 

(2) Title information (20 words). 

WRITE (LB3D) (TITLE (I) , I»l, 20) 

(3) The number of: diffusion nodes, arithmetic nodes, and total 

nodes; followed by a temperature value for each node; then 
a capacitance value for each diffusion, if any. 

WRITE(LB3D)NND,NNA,NNT, (T(l) , I=1,NNT) 

IF (NND . GT . 0) WRITE (LB3D) (C (I) , I-l, NND) 

(4) The total number of conductors followed by a conductor value 
for each one. 

WRITE (LB3D)NGT, (G(I) , I=1,NGT) 

(5) The total number of user constants are followed by the 50 control 
constant values; then the user constants values, if any. 

WRITE (LB3D)NCT, (FIXC (I) ,1=1,50) 

IF (NOT . GT . 0 ) WRITE (LBJD) (K (I) , 1=1, NOT) 

(6) The total number of arrays and the overall length of the array 
data; then the array values, if any. 

WRITE (L B3D) NAT , LENA 
IF (LENA.GT.O)WRITE(A(I) ,I=1,LENA) 
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(7) The lengths of the first and second pseudo-compute sequences, 
followed by the data for the first pseudo-compute sequence; 
then the data for the second pseudo-compute sequence, if any. 

WRITE (LB3D)LSEQ1,LSEQ2, (PI (I) , I=1,LSEQ1) 

IF (LSEQ2 . GT. 0) WRITE (LB3D) (P2 (I) , 1=1 , LSEQ2) 

Note that (3), (4), and (7) above apply only to a thermal problem. 

4.4.2 LB4P - Program F0RTRAN "Tape" 

This "tape" is especially formatted in 507 word blocks as 
required by the F0RTRAN compiler. Where: 

WORD 1 on the first block of each routine contains the name 
of the routine. 

1>DRD 2 contains the integer number of card images in the 
block. 

WORDS 3 - 506 contain the card images 

WORD 507 is set to 40 except on the last block of each 
routine where it is set to -0. 

4.4.3 INTERN - Preprocessor Scratch **Tape" 

Generally INTERN is used to pass card images to subroutine 
BLKCRD under a 14A6 format. 

4.4.4 LUTl - Dictionary "Tape" 

This "tape" contains a list of the actual SINDA numbers in a 
relative order. That is, the actual node ninnber corresponding to the 
kth relative node mnnber is the kth item of the node number dictionary. 
The format of this "tap^' is; 

(1) The total number of nodes, followed by an actual node ntaaber 

for each node. ^ 

WRITE (LUTl) NNT, (NN (I) , 1=1, NNT) 

(2) The total number of conductors, followed by the list of actual 
conductor nttmbers. 

WRITE (LUTl)NGT, (NG(I) , 1=1, NGT) 
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(3) The maaber of user constants, the total nuonber of constants, 
followed by a list of the acttml constant numbers. 


WRITE(LUT1)N0C,NCT, (NK(I),I=1,NGT) 

(4) The total number of arrays followed by a list of the actual 
array numbers, then the total number of arrays followed by a 
list of the length of each array. ^ 

WRITE(Urri)HAT, (NA(I), I*1,NAT) 
miTE(OJTl)KAT, (IA(I),I=1,NAT) 

4.4.5 LUT3 - Parameter Runs “Tape* * 

This ”tape" contains some data from the original problem. It 
is required by the initial parameters capability. The format 
of "tape" LUTS is: 

(1) The original title. 

WRITE(LUT3) (TITLE(I) , 1=1,20) 

(2) A list of original temperature and capacitance values. 

(3) 

<4) 

(5) 


WRITE(L0T3)HND, (T (I) , 1=1 , NUT) 

IF (NND.GT. 0) WRITE (LUTS) (C (I) , 1=1, NND) 

A list of the original conductor values. 

WRITE(LUT3) (G(I) , I=1,NGT) 

Lists of the original fixed and user constants. 

WRITE{LUT3)NUC,NCT, (FIXC(I) , 1=1,50) 

IF(NCT.GT. 0)WRITE (LUT3) (K(I) , 1=1, NOT) 

The original array values. 

WRITE(LUT3)NAT, LENA 

IF(LENA.GT.0)V’RITE(LUT3) (A(I) , i-l,LENA) 
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4.5 


Overlay Structure 

The SINDA preprocessor has an overlay structure composed of a 
main link (designated LINKO below) which is always in core and five sublinks 
(designated LINKl, LIIIK2, LINKS, LINK4, LINKS, and LINK6 below) which over- 
lay one another as they are brought into core. 


LINKO 


PREPK0 

BLKCED 

FINDBM 

SQUEEZ 

STFFB 

WRTBLK 


LINKl 


C0DERD 

C0NVRT 

DATARD 

ERRMES 

GENUK 

INC0RE 

N0DEDA 

QDATA 

RELACT 

SETFMT 

TYPCHK 

WRTDTA 

WRTPMT 


LINK2 


LINKS 


LINK4 


LINKS 


LINK6 

PSEUD0 

PCS2 


GENLNK 


PRESUB 

ALPINT 

MXT0FN 

SEARCH 

SINDA4 


SPLIT 

SKIP 


EDIT 










Note that the first subroutine listed above in each of the sub- 


links serves as the driver for that sublink and it is also the subroutine 
called from PREPR0. 
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READS & PR0GESSES 
DATA BL0CKS 


Another approach to overlay specification is to think of each 
link as a functional unit, hence the graph below. 
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DATA 
















4.6 


Structure of Pseudo-Compute Sequences 
Descriptions 


4.6.1 

The use and structure of the two pseMo-Compute sequences 
generated by the SINDA preprocessor appear to be rather confusing and 
mysterious. The term "pseudo" itself leads to immediate interpretation 
difficulties. Suppose that an element Gj^ between nodes i and j is to be 
identified for the ith node; by specifying explicitly i,j» and k the 
element is completely defined. On the other hand, SINDA explicitly 
specifies J and k but i is implicit in the D0-L00P; hence, the description 
pseudo-compute sequence (PCS) arises. Confusion also arises from the lack 
of information regarding the need for the (PCS) and the difficulties in 
reading the packed information. In short, the PCS as used in SINDA is 
simply two lists of relative numbers which are ordered in a specific 
manner. The two lists of relative numbers form the heart of the PCS, 
although other information pertinent to the computation must also be 
considered. 

The PCS is necessary because the data as input by the user does 
not lend itself efficiently to the computational capabilities of F0RTRAN. 

As a result, the preprocessor scans the user input data and places the 
relative numbers (F0RTRAN addresses) into an array in the order in which 
the data will be used at a later time by the user selected numerical 
solution routine. 

Packing of the data is a technique that conserves computer 
storage by placing two or more pieces of information in one computer word. 
This allows the user to execute a larger problem than the one that can 
be accommodated if the traditional one computer word for each piece of 
information approach. The penalty for this larger problem capability 
is an increase in execution time required for the extraction of information 
each time it is used. 

4.6.2 Structure of PCSl 

The first PCS, designated PCSl, contains the following information 
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bit 0 

last G for this node flag 

bit 1 

automated C option flag 

bit 2 

automated G option flag 

bit 3 

radiation G flag 

bit 4 

Q from source data flag 

bits 5-20 

relative G niamber 

bit 21 

1 way G flag 

bits 22-35 

relative adjoining node ntimber 


The 36 bits of each computer word are numbered 0 through 35 from left to 
right. All of the 1 bit flags are set such that 0 means N0 and 1 means 
YES. 

PGSl is stored in an array named NSQl and is ordered by rela- 
tive node number. That is, for relative node number 1 the conductor data 
is scanned and each time a conductor connected to node number 1 is 
encountered the PCS! information is stored in NSQl. When all of the 
conductor data has been scanned, the 0 bit of the latest word of PCS! 
information is set to 1. This process is repeated for relative node 
numbers 2, 3, etc., until all diffusion and arithmetic nodes have been 
processed. PCSl will be formed as either long or short as specified by 
the user on the BCD 3TEEEMAL card. This option is applied to diffusion 
nodes only since arithmetic nodes are always formed under the long option 
The difference between the long PCS and the short PCS is that iba the 
long PCS, each conductor will be listed twice, whereas in the short PCS 
each conductor will be listed once. This assumes the conductor connects 
two diffusion nodes. If one or both of the nodes is arithmetic, then 
the conductor will be listed twice> and if one of the nodes is a boundary 
the conductor will only be listed once. For example, given conductor 
number k which connects diffusion nodes i to j , where i < j , The long 
pcs would contain the k, j information for the processing of node i, and 
the k,i information for the processing of node j; whereas, the short 
PCS would only contain the k,j information. The short PCS thus has the 
advantage of requiring less computer storage than the long PCS, but a 
block iterative method (refer to Section 5.2.2) must be used; in general 
the short PCS requires more iterations to converge than the successive 
point iterative (refer to Section 5.2.2) method which requires the long 
PCS. 
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4.6.3 


Structure of PCS2 

The second PCS is designated PCS2. The following information 
is stored whenever bit 1, bit 2, or bit 4 of PCSl is set to one. 

bits 0-4 automated option code 

bit 5 not used 

bits 6-21 F0RTRAN address of the array or relative 

constant number. 

bit 22 not used 

bits 23-35 relative constant number 

If the automated option is a doublet type, like DIV, and therefore requires 
two words to store the information, the automated option code on the 
second word is set to zero. In the event that more than one of the 
flag bits (bits 1, 2, or 4) of PCSl is set to one, then the following 
order is imposed on PCS2: the capacitance information is stored first, 

the source block information second and finally the conductor information. 

The PCS2 information is stored in array named NSQ2. This 
array is the same under the long PCSl or the short PCSl since automated 
conductors are only flagged on their first encounter. 

4.7 Other Information 

This section contains miscellaneous information that may be of 
Interest to the user. 

4.7.1 Subroutine Lengths 

The storage required by a particular routine will vary depending 
on the t 3 rpe of computer and the system being used. The routine lengths 
given in Section 4.2 are based on compiler listings made on 23 January 1971 
at Jacobi Computation Center.* The machine is a UNIVAC 1108 with a highly 
modified system. The numbers represent the sum of the computer storage for 
computer instructions, constants, and simple variables. 

4.7.2 Maximum Thermal Problem Size and Maximum Data Value Size 

A short formula for estimating the maximum thermal problem size 
that can be run on SINDA, and a list of the maximum size of the various 
data values is given below. 

* Now called Computation and Systems Corporation, Los Angeles, California. 
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Estiaatlon of Maximum Problem Size 
HNT + 3*NGT + NCT + 4*NA0 < LENBKI 
where, 

HOT is the total liumber of nodes . 

NGT is the total number of conductors. 

NCT is the total number of constants (user 

constants plus literals from automated options) . 

NA0 is the number of automated options specified. 

- LENBKT is the length of the dynamic storage array as 

set in routine PREPR0. 


Maximum Size of Data Values 

Actual node number 

Core storage 

33 

2^-1 

Print out 

999,999 

Relative node number 
Core storage 

16,383 

Temperature 

Core storage 

± 10^^ 

Capacitance 

Core Storage 

+ 10^® 

Relative conductor number 
Core storage 

235-1 

Actual user constants number 
Core storage 

32,767 

Automated options 

16,383 

Relative user constants number 
Core storage 

32,767 

Automated options 

8,191 
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User constant values 
Integer 

Floating point 
Alphanumeric 

Actual array nxanber 


Floating point 
Alphanumeric 


±10“ 

6 eharacters 


Core storage 

2^^-l 

Automated options 

16,383 

Print out 

99,999 

Relative array number 

235-1 

Core storage 

Automated options 

65,535 

Print out 

99,999 

Array values 

± 235-1 

Integer 


±io3« 

6 characters 


Note that some of the maxima, such as the relative conductor 
35 

number of 2 -1, are strictly academic since the djnaamic storage array 

is considerably smaller than the Indicated maximum data value size. 
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5. REVIEW OF LUMPED PARAMETER EQUATIONS AND BASIC NUMERICAL SOLUTIONS 

The use of SINDA as meatioiied in a previous section is based on 
a lumped parameter representation of a physical system.^ Thus SINDA solves 
ntjmerically a set of ordinary (in gmeral nonlinear) differential equations 
that represent the transient behavior of a lumped parameter system or a 
set of nonlinear algebraic equations representing steady state conditions. 
Nimerous numerical solution techniques are reported in literature; a few 
of these are listed in the Reference Section. These numerical methods 

are based on finite difference algorithms as opposed to finite element 
methods which have received considerable attention recently. For 

problems that are generally encountered in spacecraft thermal design » use 
of the finite element method appears to be inappropriate because of the 
nonlinearity presented with radiation heat transfer and because of 
complex geometric configurations. 

Variations of the basic finite .difference algorithms are numerous 
because no single numerical solution technique is optimum for all the 
endless types of thermal problems that can be encountered. Furthermore, 
because of the nonlinearity of the problems, a specific set of criterions 
to indicate solution accuracy and stability is not. available and does not 
appear to be forthcoming. As a result, the user is placed in a rather 
awkward and confused position of not knowing which subroutine to use if a 
choice is available. Some thermal analyzer-type computer programs allow 
no choice, as a result, user decision is not necessary. SINDA represents 
a computer program at the other extreme of user decision flexibility by 
providing a number of numerical solution methods. 

The intent of this Section 5 is to review and formulate the basic 
numerical solution methods with the presentation (from an engineering stand- 
point) of the characteristics of each SINDA numerical solution routine 
deferred to Section 6. In addition to place the use of SINDA in a proper 
perspective relative to accurate temperature prediction of a physical 
system, difficulties associated with lumped-parameter representation are 
discussed here. 
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5.1 


X’troped Parameter Representation 


Reduction of a distributive (physical) system to a lumped system 
which can be represented as an equivalent thermal network is a rather impor- 
tant phase of thermal analysis. From a temperature accuracy standpoint 
lumping (or nodalization) of the physical system may be far more important 
than a numerical solution technique that is used in a computer program. The 
latter is often given undue attention with apparent ignorance of other error 
sources which may be far more important. A general discussion on lumped para- 
meter representation is not intended for presentation here since the subject 
material is extensively covered in technical literature, but it is convenient 
for continuity to Indicate basic considerations. 

For simple geometries and linear problems, it is rather 
straightforward to solve the partial differential equations of the type. 


where. 



aV^T + 


Q 


thermal diffusivity (k/C) 

temperature 

source 



(two 


dimensional) 


(5.1-1) 


Numerous analytical solutions of (5.1-1) for different types of 

3 0 31 

boundary conditions and geometries are available. * Finite difference 

al g orithms formed directly from the partial differential equations are 
also abundantly reported in literature.^^ These finite differenee 
formulations were generally developed for well-defined geometries and 
symmetrical discretization. For these problems, the so-called nodal 
connections or resistances are immediately available and, in general, 
automatically generated by the computer program. Thus, the need for a 
lumped parameter representation does not exist. For these types of 
problems, inaccuracies due to truncation and solution stability are 
specifically established. . ' ' 


For complex geometries and nonlinear problems such as those that 
include thermal radiation exchange, analytical solutions of thermal 
problems are limited, As a result, it is a common practice because 

of practical considerations to nodalize a physical system directly with- 


5 - 2 



out undue consideration of inaccuracies. Thus, a user merely represents the 
heat flow between two connecting nodes by using the basic network building 
block, 

Hi “ ^^i " (5.1-2) 

where, represents an effective resistance between adjoining nodes 

1 and j. 


It should be particularly noted here that SUnXA. employs the 
concept of conductance in lieu of resistance which is common with most 
network-type computer programs. Thus the heat flow is represented as: 


qy - GyCTi - Ij) 


(5.1-3) 


where, is the conductance between node i and node j. 


The proper value of (or R^^)for an arbitrary nodalization 
is (and should be) of concern to the user but because of the multitude of 
variables that must be considered, any discussion here would be incomplete. 
An excellent article on asymmetrical finite difference networks is 
presented in Reference 35. 


5.1.1 Some Thoughts on Ltimped Parameter Errors 

Reduction of a physical system to a topological model consisting 
of a network with resistors and capacitors requires considerable engineering 
judgment, tfore often than not, nodal size of a model is governed by budget 
and schedule constraints. As a result, the discrete areas larger than 
desired are often used. This does not necessarily mean, however, that the 
use of a large ntanber of nodes will always yield realistic results since the 
uncertainties of the input parameters can be appreciable.^^ 

Spatial truncation errors are controlled by selecting the grid size 
so that nonlinear temperature distributions lie within required accuracy by 
linear interpolation between nodal points, and that variation of temperature- 
dependent properties over the volume of each node is within required limits 
of the average values determined for the nodal point. The assximption of 
linear temperature distribution assumed for the lumped— parameter equations 
(equation 5.1-6) leads to a spatial truncation error of the order 0(Ax ) only 
if all nodes are symetrically located. If a non-uniform grid is used, the 
accuracy of computation is only 0(Z!ix). Spatial truncation errors are thus 
Inherent in the mathematical model and beyond user control once Inputted into 
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SINDA. For a two-dimensional problem with symmetrical grids, the spatial 
truncation error can be expressed for typical explicit and implicit methods 




" "2 3/“ 12 3 ^ 4 -^ 


O(Ax^) 


(5.1-4) 


Teiiq»erature distribution other than linear can also be formulated however, 
most thermal analyzer-type computer programs such as and including 

SUn}A are based on the linear assTjmption. 


Time truncation errors are directly dependent upon the time-step 
since the error for the typical er^llcit and implicit method is. 


E 


At 3^T 


3t" 


(5.1-5) 


Normally the time step is dependent upon a particular criterion chosen by 
the user. A more detailed discussion on user control of the time step 
will be given for each numerical solution routine within the SINDA sub- 
routine library .® * ** 


Another approximation error which is due to discretizing is the 
assumption of constant radiosity for the discrete areas. Inaccuracies can 
) be expected to affect the level and distribution of temperature. The 
analysis of thermal radiation exchange has received considerable attention 
in recent years because of its importance in spacecraft thermal design. '^'’* 
The Influence of non-uniform local heat flux on overall heat transfer 
between a gray differential area parallel to a gray infinite plane is 
examined in Reference 43; the assumption of uniform local heat flux appears 
to be reasonable for this geometry and for the evaluation of the overall 
heat flux calculations. A method of 2 Uialysis suitable for engineering 
applications is developed in Reference 50 for computing local radiant flux 
and local temperature of opaque surfaces in a space environment. A study 
evaluating the validity of commonly used simplified methods of radiant heat 
transfer analysis is reported in Reference 48. A study directed at 
improving the understanding and prediction of orbiting spacecraft thermal 
performance is presented in References 46 and 49. A method presented in 
Reference 51 provides a means of evaluating the uncertainties associated 
with thermal radiation exchange. For an excellent status review (as of 
1969) on radiation exchange between surfaces and in enclosures, the reader 
should consult Reference 52. 


5-4 



The above discussion merely serves to indicate that considerable 
care must be given when nodalizing a physical systas and that the numerical 
evaluation of the finite difference equations must be considered from the 
total temperature error context. This means that user attention to a given 
numerical solution must be placed in a proper perspective. 

5.1.2 Lumped-Parameter Equations 

Using the network building block as expressed by equation (5.1-3) 
the lumped parameter system is identified as a set of ordinary non-linear 
differential eqxjations by taking a heat balance as an ith node. 




P 

Z 

j=l 


a,, 

13 


P 

(T. - T.) + Z Ob 
J ^ j=l 


(T^ 

j 



(5.1-6) 


i = 1,2,...,N (number of variable temperatures 
Tj = constant, N < j j< p 

where, « the ith nodal capacity which may be a function of 

t^perature 

q^ = the heat into node i and may be a function of time and 
temperature (impressed) 

a,, = the conduction coefficient between nodes i and 1; it 

ij 

may be a function of time and temperature 
b^^ = the radiation coefficient between nodes i and j; it 

may be a function of time and temperature 
O * Stefan-Boltzmann constant 


Coefficients a^^ and b^^ are SIHDA input quantities with the 
temperature factor of equation (5.1— 6) calculated Internally by the program. 
®obh a^j and b^^ may be variables. Conductance updating is a subject for 
discussion in a later paragraph. The user requirement to input the 
coefficients, a^^ and provides considerable program flexibility, but 

at the same thue tiser generation of these input quantities presents , in 
some instances, rather difficult engineering judgment decisions. 


Radiation coefficient b.. is, in essence, a radiation inter- 

*3 5S J 

change factor, (also known as script F) between nodes i and j . 

Generation of this q\iantity analytically can be quite difficult and 
inaccurate. A number of methods and computer programs (see, for example. 
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Reference 56 ) are available for evaluating the shape factors which 
represent an Important part of detenulning script F, A direct generation 
of script F is normally through the use of the Monte Carlo technique, **^ 
but a recent development utilizes a matrix formulation for determining the 
script F in an enclosure containing surfaces with arbitrary emission and ' 
reflection characteristics.®^ An experimental technique is reported 
in Reference 39 . 


5.2 Basic Finite Ojf f erenee Formulations 

The various numerical solution techniques differ in the finite 
difference formulations for the time-derivative (refer to equation 5.1-6); 
since the thermal equation is of the parabolic type, the transient heat 
transfer problems are of the initial value type. This means that at some 
time point,- 1 “ nAt (n is the number of time steps. At) all values of T^ 
are known. Thus, 

j 

At , (5. 2^-1) 

^n,n+l 

i “ l»2t • • • 



where represents the time interval between t “ nAt and t 


(n+l)At 


It is apparent that the selection of the proper value to (dT. /dt)t 

i,n n,nTi 

cannot be explicit and its selection identifies one numerical method from 
another. Although many finite difference foirmulations of the parabolic 
differential equation are available, two general classifications are 
commonly denoted as explicit or implicit. These numerical methods are 
well-documented in literature; the reader should refer to Ref erence 12 
for a comprehensive discussion on various finite difference approximations. 
Explicit methods also discussed in References 14, 17, 19 and 20, among 
others, are step-by-step in time and equations. 


Explicit methods include: 

(1) Forward-difference explicit approximation^^* 


This is an Euler method that computes temperatures in a 
step-by-step fashion. The requirement of stability 
places an upper limit on the time increment. SINDA 
subroutines CMFRUD, CNFRDL and CNFAST fall within this 
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category. CNFAST is a modified versioa of CEFRW which 
allows the user to specify the minimum time step to be 
taken. Refer to Sections 6.3.1 and 6.3.2 for details. 

(2) Dufort-Frankel approximation^ » 

The Dufort and Frankel finite difference formulation is a 
three level formula that appears to be unconditionally 
stable. SIKDA subroutine CNpUFU uses the Dufort-Frankel 
finite difference algorithm (refer to Section 6,3.4 for a 
detailed discussion). 

(3) Exponential approximation ^ 

The exponential approximation is found by integrating the 
heat balance equations after making linear and constant 
coefficient assumptions. This method is unconditionally 
stable for linear systems but may be unstable for some 
types of nonlinear problems. SlNDA subroutine CNEXPN 
employs this method and is discussed at length in Section 6.3.3., 

(4) Alternating direction approximations^’ 

This technique employs two formulations, one on odd time 
levels and the other on even time levels and is uncondi- 
tionally stable. 

The implicit finite difference formulations require a simultaneous 
computational procedure, In addition to Ref erence 12 , Implicit methods 
are also discussed in References 8, 10, 17, and 20, among others. Implicit 
methods include: 

(1) Backward difference implicit approximation^^ 

The backward difference weights only the flux terms at 
t = (n+l)At. As a result, the method is stable for all 
V£ilues of At. SINDA subroutine CNBACK employs this method 
and is detailed in Section 6.4.1. 

(2) Crank-Nicolson approximation® 

The Crank-Nicolson method uses the arithmetic average of 
the heat flux at the two time levels, t = nAt and 
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t * (n+l)At. The methcd Is unconditioBally stable. SINDA 
CN7WBK uses this method and is discussed in Section 6.4.2. 

Steady state analysis also requires an implicit method of 
solution. SINBA steady state subroutines are called CIHDSS, ClhBSL and 
CINDSM which are detailed in Sections 6.5.1, 6.5.2, and 6.5.3, 

5.2.1 Forward Finite Difference Explicit Method 


By replacing the first derivative of temperature with respect to 
time, dT/dt, with the forward first difference quotient, equation (5,1-6) 
becomes. 


C ^i.n+1 “ ^i.n^ - 

At 'll 


P 

Z 

j-1 




T- „) 

x,n 


P 

Z 

j-1 






t ■ nAt 

1 « 1,2, • ,N 


j.n 


constant, N < j < p 


where, the second subscript on T represents the time level such that 


Equation (5.1-6) is represented in the form expressed by equation 
(5.1-3) by letting, 

(T^ + Tj)(T^ + T^) (5.2-2) 


It is interesting to note that the finite difference form of (5.1-6) (and 
thus 5.2-D represents a second central-difference quotient of V T (refer to 
(5.1-1). 

The computational procedure for the forward difference formulation 
is rather straightforward since only a single unfcaown temperature at each 
time Step, T = nAt for each equation is present. Note that the averaging 
of dT^/dt) assigns a weighting factor to the heat flux terms only (terms on 
the right side of equation (5.1-6) at t = nAt). Along with the computational 
simplicity, however, is the stability constraint which places an upper 
limit on the time increment. At, that can be used in the numerical pro- 
cedure. The stability criterion for the explicit finite difference method 
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is (for the Eost limiting node) , 


P 

At < C./ 2 G., 


(5.2-3) 


A modified stability criterion that allows for a larger time step 
which results in a conditionally stable temperature for the most limiting 
node is reported in Reference 23. Since the stability criterion will 
govern the maximum time step that can be used, it is thus particularly 
Important that a user gives some attention to those factors that compose 
the condition of stability when nodalizing a physical system. 


■ In the discussion presented so far, arithmetic nodes (nodes with 
no heat capacity) have not been mentioned. Normally, the computational 
procedure treats arithmetic nodes separately from the diffusion nodesj 
arithmetic-node temperatures are solved Implicitly. Detailed discussion on 
the general procedure will be presented in a later paragraph as well as in 
Section 6 which discusses the various SINDA numerical solution routines i 

5.2.2 Implicit Finite Difference Method 


The implicit difference equations can be constructed for heat 
transfer problems in many ways (see, for example. References 12 and 20). 

Replacement of equation (5.1-6) with the backward time difference 

yields. 


(T. . - - T. ) 

i,n"bl i ,n 

^i At 


P 


jfi ^ij " '^i,n+l 


) + I ah,. (T^ 

ij j »n+l i,n+l 


(5.2-4) 


i * 1,2,...,N 
Tj * constant , N < j _< p 

' ■'i 


The computational procedure for the backward difference formula- 
tion must necessarily be re— iterative because of the need to solve a set of 
simultaneous non-linear equations. 


In view of the importance of iteration techniques (such as method of 
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successive approximation), it may be of interest to formulate equation (5.2-4) 
into an interative form. If we let use equation (5.2-2) in equa- 
tion (5.2-4) and solve the resultant expression for n+1* yields the 

recurrent equation for a given time increment. At, and time-step, n. 


i,lcH 


1 , T . , + £ G., , T, , + q. , 

i,k i,h X3,k 3,k ^x,k 


'^i.k + jfj ®«.k 


i = 1,2, ...,N 


(5.2-5) 


where, = C^^^^/At 




,^ij,k °'^ij,k ^'^j,k ’^i,k^^’^j,k ^i,k^ 


(5.2-6) 


T , » constant, N < j ^ P 

J 

k “ kth iteration (note that C 


i,k» '^i,k» “ij •'ij 
to be updated every iteration; SINDA routines update these 

quantities once each time-step) 


a . , and b . . are shown 


The iterative pattern is initiated by assuming "old" temperatures 
(Xi j, and ^ on the right side of equation (5.2-5) to evaluate a "new" 
set of temperatures (T^ on the left side of the equation (5.2-5); 

this single set of calculations represents an iteration. By replacing all 
of the "old" temperatures (T^ ^ on the right side of eqxaation (5.2-5) 
with the just calculated "new" set of temperatures (T, ,^, ), a second 
iteration can be made. The iteration procedure is continued until a 
termination criterion such as the number of iterations or the maximum 
absolute difference between T^ and T^ is less than some prespecified 
value has been satisfied. It should be noted that 
shown to be updated every iteration. This iterative process is termed 
"block" iteration since the "old" temperatures on the right side are 
replaced in a "block" (a set of temperatures) fashion with the "new" 
temperatures . 


Another iterative technique is to utilize on the right side of 
equation (5.2-5) each "new" temperature as soon as it is calculated. This 
iterative method is termed "successive point" iteration and appears to • 
yield solutions about 25% faster than the "block" iteration method. 

Equation (5.2-5) can be expressed in a "successive point" form 
as follows: , 
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where. 


%k "i,k ®ij,k ■^j,k+l ®ij,k ^j,k ^ ^ 

^i,k+l ~ _ P 

^ ®ii k 

1 j=l 

=«.k - “ij.k + '^13 ,k + n.k><^3.i ■" ^i.k^ 


^ (5,2 


{it = k if j > i and it » fcfl if j < i) 

Tj ^ “ constant^ N < j £ p 

k “ kth iteration (note that C. ^ i.* * 44 * ^44 *^® shown to be 

X,K X,iC Xj X3 

updated every iteration) 


. The iterative method as used in SINDA follows a fixed, pre- 
determined sequence of operations in contrast with a relaxation procedure 
which is also one of successive approximations but is not processed out in 
a predetermined sequence. The relaxation procedure seeks and operates on 
the node with the maximum temperature difference between the "old" and 
the "new," Prom a programming standpoint, the search operation requires 
as much computational time as the temperature calculation itself. 

5.2.3 Steady State Method 

Standard steady state equations follows directly from equation 
( 5 . 2 - 5 ) for block iteration or from equation (5.2-7) for successive point 
iteration by letting = 0 in these equations. The comments made in 
Section 5.2.2 are equally applicable here. 

5.2.4 Some Comments 

The finite difference expressions presented in this Section 5.2 
represent standard formulations and thus do not show computational tech- 
niques and artifices which are used, some more or some less, in all programs. 
SINDA ntoaerical solution routines contain many computational features (many 
original with J. D. Gaski) which enhance problem solution. The various com- 
putational aspects of the ntnnerical solution methods as used in the SINDA 
routines are discussed in rather lengthy detail in Section 6. 
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6. SINDA NmiERICAL SOLUTION ROUTINES 

Objective and Presentation Format 

SINDA has available to the user a number of numerical solution 
routines which ^iploy various numerical methods. A brief description of 
these routines with the required SINDA input quantities and format are 
contained in the SINDA users manual;^*** a general review of numerical methods 
was presented in Section 5, Unfortunately, the brief description is not 
sufficient for a casual SINDA user to make a selection decision from among 
several routines that are available and for a serious user to fully under- 
stand the computational procedure as well as to xinderstand the role of the 
various control constants that are employed in each routine. 

It is the intent of this section to fill wherever possible and 
practical the description void that presently exists with the numerical 
solution routines by detailing the characterisitcs of each. It is not the 
intent here to provide sufficient detailed Information for a user to make 
modifications and/or additions to the existing subprograms, but rather to 
provide information that will aid the user in assessing the various 
numerical solution routines and in evaluating the numerical results. 

Each of the numerical solution routines is detailed from a x 

theoretical as well as from a computational standpoint. Control constants 
and their role are described and Indicated in a step-by-step verbal flow 
computational procedure. Details of many of the numerous computational 
checks have purposely been omitted because of the complex interactions. 

Minute details of each routine can be obtained only from the individual 
computer listings; a computer listing of each of the SINDA nvnnerical 
solution routines is presented in Appendices A, B and C. General computa- 
tional procedure and features that apply to most, if not all of the SINDA 
numerical solution routines are assembled in a single section (6.2) in 
order to eliminate undue repetition. The description of each routine is 
heavily dependent upon, and coupled to, the general description of 
Section 6.2. The routines have been categorized as steady state or 
transient with- the latter subcategorized as explicit oi: implicit in order 
to allow for an orderly presentation as well as to simplify future additions. 


6-1 



6.2 


General Computational Procedure and Features 

Ea<ch of the SINDA numerical solution routines employs a particular 
finite difference approximation of the liunped parameter heat balance equa- 
tions. In spite of the uniqueness of each routine, portions of the computa- 
tional procedure used in each are similar. Also, many of the routines have 
identical features such as the acceleration of convergence and the use of 
control constants. As a result, it is convenient to place in this section 
repetitious material. In some instances material presented here is repeated 
in the discussion of a particular numerical solution routine. 


6»2.1 Order of Computation 


It was reported in Section 3.5 that the order of computation depends 
on the sequence of subroutine calls placed in the EXECUTI0N block by the 
program user. No other operations block is performed unless called upon 
by the user either directly by name or indirectly from subroutines which 
internally call upon them. Nimierical solution subroutines internally call 
upon operations blocks VARIABLES 1, VARIABLES 2, and 0UTPUT CALLS. The 
internal order of computation for these numerical solution routines is 
similar with the primary difference between one routine and another being 
the finite difference approximation employed in a particular routine, A 
flow diagram indicating the general order of computation for the numerical 
solution routines is depicted in Figure 6.2-1. 


6. 2. 1.1 Finite Difference Algorithm 

Although each of the SINDA numerical solution routines employs a 
particular finite difference approximation which is detailed for each 
numerical solution routine, the computational pattern is similar. Within 
the box depicted as SFDA | in Figure 6.2-1, solution of the finite 
difference algorithm occurs. The computational sequence for transient 
solutions follows one of two patterns: (1) one for explicit finite 

difference methods; and (2) one for implicit finite difference methods; 
steady state solutions follow closely the implicit pattern. Both numerical 
flow pictures are depicted in Figure 6.2-2; details within the flow pictures 
are different for each routine and are described separately under the indi- 
vidual SINDA numerical solution routines (refer to Sections 6.3 - 6.5). 
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OPERATION 
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j GTS 

] 

Calculate time step 

1 VARBLl 

] 

Variables 1 operation 

1 SFDA 

3 

Solve finite difference 



algorithm 

1 VARBL2 1 

Variables 2 operations 

1 0UTCAL 

3 

Output calls operations 

1 MTC 

3 

Modify time control 

t El 

3 

Erase Iteration 

Check 


Reverse direction if 



Backup nonzero 



Relaxation criteria not 
not met 



Time or temp change 
too large 



Backup nonzero 

Y 


Not time to print 

§ 


Problem stop time 
not reached 

A 


Implicit routines, 
iteration loops 

B 


Implicit routines, 
once per time-step loop 


Figure (6.2-1) General Order of Computation 
for Numerical Solution Routines 
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For CINDA -3G users, it should be noted that the updating of properties occurs 
within the numerical solution routine after VARIABLES 1 call. CINDA-3G 
evaluates the variable properties before VARIABLES 1 call. 


Figure 6.2-2. 


Numerical Computational Pattern for Explicit 
and Finite Difference Algorithms 
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6. 2. 1.2 Updating of Optionally Specified Properties 

Optionally specified properties are defined here as those items 
which result in pointers being set in the second pseudo compute sequence 
(refer to Section 3,3.4). The term optional refers to mnemonic options that 
are available for different types of variable properties.^* ^ The properties 
are updated in all SINDA n'umerical solution routines the same way. This definition 
is used here in lieu of stating that optionally specified properties are time 
and/or temperature varying properties since source date may be specified to 
be constant. The pointers are set by one or more of the following user input 
quantities; 

(1) All capacitances t C^, specified as f (T) or f (t,T) in NODE DATA 
“BLOCK: 

(2) All data, q^, entered in the SOURCE DATA BLOCK: 

(3) All coefficients, Gj^, specified as f(T) or f(t,T) in CONDUCTOR 
DATA BLOCK. It should be noted here that the term coefficient as 
used here requires amplification. The conductance, , may be 
for conduction or for radiation; that is. 


\ 


a^. (for conduction conduetance) 


®ij " ®^ij ^^i + radiation conductance) 

- Gy. (T^ + Tj )(X^ + Tj) 

Thus, note that the calculated conduction conductance is 
Identical to the updated 6^^, whereas for the calculated radiation 
conductance only ab^^ is equivalent to the updated Gj^. 

The type of optional properties is identified by the integer stored 
in the first six bits of the second pseudo compute sequence which indicates 
to the program which option is in effect. Optional property types are listed 
and described for the three categories of input quantities in Table 6.2-1 
for capacitance. Table 6.2-2 for Impressed source, and Table 6.2-3 for 
coefficients with the definition of symbols listed in Table 6.2-4. 


6.2.2 Operations Blocks 

In a previous paragraph, it was mentioned that the sequence of 
subroutine calls placed in the EXECUTION block by the user determines the 
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TABLE 6.2-1 


OPTIONALLY SPECIFIED CAPACITANCE EXPRESSIONS 


Option 

l2Ee 

Expression 

SIV 

1 

= F(A^:T^) 

DIV 

2 

C^ = FKaJsT^) + F2(A2:Tp 

DIV 

3 

* F1(L) + F2(A":T^) 

DIV 

4 

* F1(A^:T^) + F2(L) 

SPV 

5 

C^ » F(A^:Tj^) 

DPV 

6 

C. * Fl(Af:T.) + F2 (a5:T.) 

X lx Z X 

DPV 

7 

C^ * F1(L) + F2(aP;T^) 

DPV 

8 

C^ = Fl(AP:Tp + F2(L) 

BIV 

9 

<=1 - F(a’>:T^. 


Notation: Refer to Table 6.2-4. 


TABLE 6.2-2 OPTIONALLY SPECIFIED IMPRESSED SOURCE EXPRESSIONS 


Option Type 


Expression 


blank 

SIV 

SIT 

DIT 

DIT 

DIT 

DTV 

DTV 

DTV 

Notation: 


1 + F 

2 “ *^i 

3 V 

4 q. = q. + Fl(Aj:t^) +F2(4:tJ 

5 q^ = q. + F1(L) + F2(A^:t^) 

6 q » q. + Fl(A^:t J + F2(L) 

X m 

7 q^^ » q^ + Fl(Aj:t^) + F2(A2:T^) 

8 q^ = q^ + F1(L) + F2(A^:T^) 

9 q. = q^ + Fl(A^:tj^) + F2(L) 

Refer to Table 6.2-4. 
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Table 6,2-3. Optionally Specified Coefficient Eispressions 
for Conduction and Radiation 

Type Egpressioa 


SIV 


1 

G. * F(A^; T ) 
k n 

SIV 


2 

Gj. = F(A^:T^) 

DIV 

(conduction) 

3 

= 1.0/[1.0/F1(A^:T^) + 1.0/F2(A2:T^)] 


(radiation) 


G^ - [FI (aJ : T^) ] [F2 (aJ : T^ ) 1 

DIV 

(conduction) 

4 

Gy. = 1.0/[1.0/F1(L) + 1.0/F2(A^:Tj)3 


(radiation) 


Gj^ = [F1(L)3[F2(A^-T^)3 

DIV 

(conduction) 

5 

Gj^ » 1.0/I1.0/F1(A^:T^) + 1.0/F2(L)] 


(radiation) 


Gj.= [F1(A^:T^)][F2(L)3 

SPV 


6 

G, » F(A^:T ) 
k m 

SPV 


7 

Gy. = F(A^:T^) 

DPV 

(conduction) 

8 

Gj, » 1.0/[1.0/F1 (aP:T^) + 1.0/F2(a|:T^)3 


(radiation) 


Gy. ^ [FKaJcT^) 3 [F2 (aP ;T^ ) 3 

DPV 

(conduction) 

9 

Gj^ * 1 . 0/ [1 . 0/Fl (L) +1. 0/F2 (A^ : T^ ) 3 


(radiation) 


\ = [F1(1)3EF2(aP;Tj)3 

DPV 

(conduction) 

10 

Gy. = 1.0/I1.0/Fl(AP:Tp + 1.0/F2(L)3 


(radiation) 


Gj^= [Fl(AP:T^)3[F2a)3 

BIV 


11 

G ■ = F(A^:T , t ) 
k ■ m* m 

SIV 


12 

Gy^ = F(A^T^) 

SPV 


13 

G, = F(A^:T.) 
k j 


Notation: Refer to Table 6.2-4; note ab^. (for radiation) 


= a. . (for conduction) 
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Table 6.2-4. Definition of Symbols for Tables 6.2-1 — 6.2-3 


Symbols 
F, FI, F2 

Gk(-aby) 


At 


m 


m 


IE 

(A^:T^) 


Definition 

Capacitance of itb node. 

Multiplying factors, either user constants or literal 
Conduction coefficient. 

Badlation coefficient. 

A literal multiplying factor. 

Heat load into the ith node, (impressed) 

Time-step 

Mean time, (TIME0 + TIMEN)/2.0 
Mean temperature, (T^ + Tj)/2.0 

Interpolated value of array A using t^ as the independent variable. 
Interpolated value of array A using as the independent variable. 


'\ b 

) (A :T., t ) Interpolated value of the bivariate array A using T. and t as 
independent variables. x m 

(A^:T^, t^) Interpolated value of the bivariate array A using and t^ as 
independent variables. 


m m 


m 


m 


Maemonic Options 


BIV Mvarlate foterpolation Variable 

DIT j^uble ^terpolation with Time as variable 

DIV Double foterpolation Variable 

DPV Double ^lynoraial Variable 

DTV ^uble interpolation with Mme and T^perature as Variables 

SIT jingle jbterpolatlon with Kme as variable 

SIV Jingle ^terpolation Variable 

SPV Mngle ^lynondal friable 

Subscripts 

1 Indicates the ith node. 

J Indicates the jth node. 

2 Indicates two (array). 
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order of computation. Operations blocks number four, EXECUTION, VARIABLES 1, 
VARIABLES 2, and OUTPUT CALLS, These operations blocks are described in 
the SINDA Users Manual^ * but their role insofar as the numerical solution 
routines are concerned may be of particular interest. 

6.2. 2.1 EXECUTION Operations Block 

The EXECUTION operations block provides the user considerable 
flexibility in the use of SINDA calls and F0RTRAN operations. Combinations 
of SINDA calls and P0RTRAN operations are innumerable since the user is 
actually programming. Now all Instructions contained in the VARIABLES 1, 
VARIABLES 2, and 0UTPUT CALLS are performed each iteration or on the output 
call interval. Thus, if an operation being performed in VARIABLES 1 
utilizes and generates non-changing constants, the operation should be 
placed in the EXECUTION block (prior to the numerical solution call) so 
that it will be performed only once and thus elitainate repetitious non- 
changing calculations. Operations of this type are conveniently performed 
in the EXECUTION operations block. Note, however that a constant impressed 
source should be placed in the optional source data block for SINDA and 
VARIABLES 1 block for CINDA-3G, 

6. 2.2. 2 VARIABLES 1 Operations Block 

The VARIABLES 1 operations block provides the user with a 
means of specifying at a point in the computational sequence, as shown in 
Figure 6.2-1, the evaluation of nonlinear network elements, coefficients 
and boundary values not considered by the various mnemonic codes utilized 
for node, conductor and source data. It is seen from Figure 6.2-1 that 
VARIABLES 1 operations occur just prior to entering the numerical solution 
phase in order to define the network completely. 

6.2 .2.3 VARIABLES 2 Operations Block 

VARIABLES 2 operations are post-solution operations in contrast 
to the VARIABLES 1 operations which are pre-solution operations as shown 
in Figure 6.2-1. VARIABLES 2 provides the user with a means to examine 
the characteristics of the numerical solution and make corrections. For 
example, the heat flow from one node to another can be evaluated or a 
temperature (s) determined without material phase change can be corrected 
to account for the phase change by using the VARIABLES 2 operations block. 
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e.2,tA OUTPUT CALL Operations Block 

The OUTPUT CALL operations block provides the user with a means of 
calling any desired subroutine with the operation performed on the output 
interval. In addition to various subroutines for printing output, several 
plotting subroutines are available.® * ** 

6.2.3 Control Constants 

Control constants number forty-nine and have alphanumeric names. 
Control constant values are communicated through program common to specific 
subroutines which require them. Whenever possible, control constant values 
not specified are set internally to acceptable values. If a required con-^ 
trol constant value is not specified, an appropriate error message is printed 
and the program terminated. Each of the SINPA niimerical solution routines 
employs a number of control constants which fall under the categories as: 

(1) user specified; (2) optionally user specified; (3) internally set by 
program; and (4) dtnnmy. These control constants are listed alphabetically 
with a brief description of each in Section 6. 2. 3.1 followed by a detailed 
description of user specified control constants in Section 6. 2.3.2; nominal 
values cf these control constants that must be specified or are optionally 
specified for each SINDA numerical solution routine are indicated in 
Table 6.2-5. Specification of these control constants is detailed under the 
discussion of each SINDA numerical solution routine. 

6.2. 3.1 Alphabetical Listing and Brief Description of Control Constants 

ARLXCA (control constant 19) 

Maximum arithmetic node relaxation temperature change allowed 
between iterations; this check occurs after each iteration. 
Specification is required for the implicit and steady state 
routines (except CINDSM) and if not specified an error 
message is printed if the number of arithmetic nodes is greater 
than zero. Specification is not required for explicit routines 
and if not specified, AELXCA is set to l.E+8. 

ARLXCC (control constant 30) 

Maximum arithmetic node relaxation temperature change calculated 
by program; ARLXCC ^ ARLXCA check is made. 
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Notation: * Indicates '’run” termination with error ffcasage printout if not specified. 



ATMPCA 


ATMPCC 


BACKUP 


BALENG 


CSGFAC 


CSGMAX 


CSGMIN 


(control constant ll) 

Maximum arithmetic node temperature change allowed between 
time steps for transient routines; the check occurs after the 
specified number of iterations. If not specified or if specified 
to be 0.0, ATMPCA is set to l.E+8. 

(control constant 15) 

Maximum arithmetic temperature change calculated by program; 
AEMPCC ATMPCA check is made. 

(control constant 12) 

Backup switch that is checked after VARiABLES 1 and VARIABLES 2 
calls. Initialized at zero. If specified to be non-zero, the 
completed time step is erased and repeated, 

(control constant 33) 

A user specified system energy balance to be maintained; this 
control constant is presently used only in CINDSM. If not 
specified, an error message will be printed. 

(control constant 4) 

Time step factor for explicit routines except CNFAST. If not 
specified or if specified to be less than 1.0, CSGFAC is set 
internally to 1.0. 

(control constant 23) 

Maximum value of C^/E ! this value aids in the checkout of 
the thermal network and is calculated only by the output sub- 
routines, CSGKMP and RCDUMP. 

(control constant 17) 

Minimum value of C./E G . . ; this value is used to limit the 

X ij 

computational time step for explicit methods of solution. 

If CSGMIN is calculated to be £ 0.0, an error message is printed. 
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CSGEAL 


DAM>A 


DAMPD 


DRLXCA 


DRLXCG 


DTIMEH 


DTIMEI 


DTIMEL 


(control constant 24) 

Allowable range between CSGMIN and CSGMAX: this control constant 

is not presently used but is included for future considerations. 

(control constant 9) 

Aritbaietic-node damping factor for all numerical solution routines; 
if not specified, or if specified to be < 0.0, DAMPA is set to 1.0. 
(Refer to equation 6.2-6.) 

(control constant 10) 

Diffusion node damping factor for implicit and steady state routines; 
if not specified or is specified to be ;< 0.0, DAMPD is set to 1.0. 
(Refer to equation 6.2-20.) 

(control constant 26) 

Maximum diffusion node relaxation temperature change allowed between 
iterations for implicit and steady state routines; this check occurs 
after each iteration. If not specified an error message will be 
printed when the number of diffusion nodes is greater than zero. 

(control constant 27) 

Maximum diffusion node relaxation temperature change calculated by 
the program; DRLXCA _< DRLXCC check is made. 

(control constant 8) 

Maximum time step allowed; applies to transient routines. If 
not specified or if specified to be ^ 0.0, DTIMEH is set to 
1-.E-8. 

(control constant 22) 

Specified time step for implicit solutions; if not specified, 
an error message will be printed and the "run" terminated. 

(control constant 21) 

Hinlm\M time step allowed for explicit routines. If not specified 
for CNPAST, an error message will be printed and the "run" 
terminated. If DTIMEO is less than DTIMEL the routines will 
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DTIMEU 


DTMPCA 


DTMPCC 


ENGBAL 


LAXFAC 


LINECT 


L00PCT 


N0C0PY 


NL00P 


tenainate with an error message, except for CNFAST which will 
do a steady state solution on the offending node. For all 
routines DTIMEL is initially set at 0.0 internally. 

(control constant 2) 

Contains time step used in computational procedure. 

(control constant 6) 

Maximiom diffusion node temperature change allowed between 
time steps for transient routines. If not specified or if 
specified to be 0.0, DTMPCA is set to l.E+8. 

(control constant 15) 

Maximum diffusion node temperature change calculated by program; 
DTMPCA £ DTMPCC check is made. 

(control constant 32) 

Calculated energy balance of the system; presently used only 
in CINDSM. 

(control constant 49) 

Specified nximber of iterations to be performed on a linearized 
system with no updating of elements during a set of LAXFAC itera- 
tions for CINDSM only; if not specified, an error message is 
printed and the "rtin" terminated. 

.(control constant 28) 

A line counter location for program output (integer) . 

(control constant 20) 

Contains number of Iterations performed (integer) . 

(control constant 34) 

Contains the no copy switch for matrix users. 

(control constant 5 ) 

Ntjmber of specified iteration loops. Must be specified for the 
steady state and implicit routines; if not specified, an 
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error message is priated and the "run” is termiaated. Optional 
specification for solution of the arithmetic nodes in the 
explicit routines; if not specified, NL00P is set to integer 1. 

0PEITR (control constant 7) 

Output each iteration if 0PEITR is specified to be non-zero; if 
not specified, 0PEITR is set at zero. May be switched on and 
off during a run. 

0UTPUT (control constant 18) 

Time interval for actiyating 0TJTPUT CALLS of transient routines; 
if not specified, error message is printed and the "run" terminated. 
May be addressed by user and modified during a run in VARIABLES 2 . 
Can be used in steady state routines for a series of steady state 
solutions. 


PAGECT 


TIMEM 


TIMEN 


TIMEND 


(control constant 29) 

A page counter location for program output (integer). 


(control constant 14) 

Mean time for a computation interval; TIMEM 
(control constant 1) 


TIME0 + TIMEN 

2.0 


New time at the end of the computational interval. 
TIMEN S' TIME0 + DTIMEU. 


(control constant 3) 


Problem stop time for transient analysis. Must be > TIME0 for 
all routines; if not, an error message is printed and "run" 
terminated. May be addressed by the user and modified during 
a run. 


TIME0 (control constant 13) 

Old time, at the start of the computational interval. Also 
used as the problem start time and may be negative; if not 
specified, TIME0 is set at zero. 
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ITEST, JTEST, KTEST, LTEST, MTEST (control constants 39, 40, 41, 42 
and 43, respectively) 

Contain dunmy integer constants. 

RTEST, STEST, TTEST, UTEST, VTEST (control constants 44, 45, 46, 47, 
and 48, respectively) 

Contain dxffinay floating point constants. 

(Control constant 31) 

Problem type indicator, 0 = THEEMAL SPCS, 1 * THESMAL LPCS, 

2 f GENERAL. 

(Control constant 35) 

Contains relative node number of CSGMIN. 

(Control constant 36) 

Contains relative node number of DIMPCC. 

(Control constant 37) 

Contains relative node ntxmber of ARLXCG. 

(Control constant 38) 

Contains relative node number of ATMPCC. 

6. 2.3. 2 User Specified and Optionally User Specified Control Constants 

The availablity of control constants which must be specified or 
which can optionally be specified provides the user with considerable 
flexibility to alter the computational criteria and hence the calculated 
temperatures. On the other hand, this flexibility presents the user with 
the problem of imputting control constant values if the nominal values are 
not suitable. An attempt will be made here to provide some guidelines 
on control constant values based on rather limited data presently available, 
but it should be recognized that suitable values to be used are dependent 
on the problem to be solved and often a trade-off must be made between 
accuracy and computational time. This normally can be obtained only 
through the use of the numerical solution routines. 
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ARLXCA (Allowable Arithmetic Node Relaxation -Temperature Change) 

This control constant must be specified for the implicit routines 
if any arithmetic node is present and for the steady state routines except 
CIHD^. For the explicit solution routines, AELXCA may be optionally 
specified; if not specified AELXGA is set to l.E+8. ARLXCA represents 
a maximum temperature change convergence criterion for the arithmetic nodes; 
ARLXCA is checked each iterative step. It is used in conjmction with 
Control constant NLO0?. Satisfaction of either NL00P or ARLXCA during any 
iterative step terminates the arithmetic node temperatures calculation for 
that time-step with computation proceeding on to the next one. Typically, 
an ARLXCA value is 0.01, but its value is dependent upon the magnitude of 
expected -temperatures. The 0.01 value tries for 5th digit accuracy for 
temperatures in the hundreds. An ARLXCA value of 0.0001 would try for 
seventh digit accuracy. Since the computer will not yield 8 digit accuracy, 
an ARLXCA value < .0001 will always result in NL00P iterations being 
performed. 

ATMPCA (Allowable Arithmetic Rode Temperature Change) 

This control constant may be optionally specified by the user for 
the implicit routines and for the explicit routines except CNFAST. If not 
specified, ATMPCA is internally set at l.E+8. ATMPCA represents an 
allowable arithmetic-node temperature change criterion between one time- 
step and another with the calculated temperature change stored in control 
constant ATMPCC. If the maximum arithmetic-node temperature change is 
greater than ATMPCA, the time-step. At, is shortened to. 

At * .95 *At (ATMPCA/ ATMPCC) 

and the arithmetic-node and diffusion-node temperatures re-set to former 
values. The computational procedure is repeated with the smaller time-step. 
Specification of ATMPCA prevents a rapid temperature change between time- 
steps with the value to be specified dependent upon the problem. Thus, the 
user should estimate the number of time-steps and the range of the tempera- 
ture to arrive at a reasonable value. For typical spacecraft-type thermal 
problems an ATMPCA of about 10“F is typical. 


6 - 17 



BACKUP (Backup Switch) 

} Control constant BACKUP provides the SINDA user with the means 

to utilize any thermal numerical solution subroutine as a predictor program. 
All of the nwerical solution subroutines set control constant BACKUP 

to zeroy just prior to the call on VARIABLES 2. Then immediately after the 
return from VARIABLES 2. a nonzero check on BACKUP is made. If BACKUP is 
nonzero, all temperature calculations for the jiist completed time-step 
are eliminated, the old temperatures (temperatures calculated at the 
previous time-step) are placed in the temperature locations and the 
control is routed to the start of the computational sequence. 

It should be noted that the user must provide the necessary 
check and criterion in VARIABLES 2 if the iteration is to be repeated. 

Thus, if the iteration is to be repeated, BACKUP must be nonzero and a 
criterion that can be met in subsequent passes established. For example, the 
criterion may require the correction of a parameter used by the network 
solution. Further, if other calls in VARIABLES 2 are not to be performed 
F0RTRAN instructions must be generated to bypass these calls. 

) It should be noted that BACKUP is sometimes checked after 

VARIABLES 1. Ifowever, for the present this use should be Ignored since 
BACKUP check after VARIABLES 1 is planned for future additions of special 
boimdary calculation subroutines. 

BALENG (User Specified System Riergy Balance ) 

This control constant is presently used in the steady state 
routine CINDSM but not in the other SINDA numerical solution routines. 

BALENG must be specified, otherwise the "run" is terminated with an 
error message printout; the value of BALENG is a criterion that represents 
an acceptable net energy balance (energy in minus energy out) of the system 
in the calculation of steady state tonperatures . A value for BALENG 
depends upon the magnitude of system energy under consideration. As a 
guideline 1/2% of the total energy into the system (including heat flow 
from the boundary) is a reasonable value. 

CSGFAC (Time Step 'Factor) 

This control constant may be optionally specified by the user 
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for the explicit routines except CHFAST and it provides the user with some 
control on the compute time-step as indicated in Section 6.2.4. If CSGFAC 
is not specified or is specified to be less than one by the user, it is 
internally set at 1.0. For subroutines CNFRWD and CNiTOL which are con- 
ditionally stable CSGFAC is a divisor; a value of CSGFAC greater than one 
is used to obtain higher accuracy. For subroutines CNEXPN, CNDUFR and 
CNQUIK, which are unconditionally stable, CSGFAC is a multiplier (refer to 
page 6-24); a value of CSGFAC greater than one is used to deerease the 
computational time. A question may be raised, why a value of CSGFAC less 
than one is not allowed for CNEXPN, CNDUFR and CNQTJIK? The reason for 
this is that it is more accurate to use CNFRWD (or CNFRDL) if a smaller 
time-step than the one associated with CSGFAC equal to one is desired. 

DAMPA (Damping Factor for Arithmetic Nodes) 

This control constant may be optionally specified for all of the 
SINDA numerical solution routines; if not specified, of if specified to be 
0.0, DAMPA is set to 1.0. In the development of the finite difference 
expressions as reported in technical literature, little (if any) mention is 
made about the so-called damping factor. The damping factor does nothing 
more than to allow a certain fraction (1.0 - DAMPA) of the "old" temperature 
(temperature at the previous time-step or iteration) to be included as part 
of the temperature change for the current time-step or iteration. The 
value to be used is dependent upon the problem and to some extent upon the 
routine. Typically, a value of 0.6 is used but a value as small as 0.01 
has been used with CIMDSL for a thermal radiation-dominated problem. In 
general, a choice for DAMPA becomes a trial and error procedure. DAMPA is 
used only with arithmetic nodes (refer to equation 6.2-6). 

DAMPD (Diffusion Node Damping Factor) 

This control constant may be optionally specified for the implicit 
and steady state routines; if not specified or if specified to be ^ 0.0, 
DAMPD is set to 1.0, DAMPD serves the same purpose for the diffusion nodes 
as DAMPA provides for the arithmetic nodes (refer to equation 6.2-21). 

DRLXCA (Allowable Diffusion-Node Relaxation Temperature Change) 

This control constant must be specified for the implicit routines 
and for the steady state routines except CINDSM. DRLXCA serves the same 
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purpose for the dlffusioii'-iaocies as control constant ASLXGA. does for the 
crithmetic nodes. Thus, the discussion on ABI^CA equally holds true for 
DKLXCA. It may be asked, why AKLXCA and DRIXCA? The reason for this is 
that it provides greater computational flexibility^ 

DTBffiH (Maximum Time-Step Allowed) 

This control constant may be opti<Maally specified for the explicit 
and the Implicit routines. DTIMEH represents the maximum time-step allowed 
during the computational process. One use of DTIMEH is the prevention of a 
single large and a single small computational time-step during an output interval 
by specifying DTIMEH as a fraction of the output interval. If DTIMEH is not 
specified, DTIMEH is set to l.OE+8. 

DTIMEI (Specified Time-Step for Implicit Routines) 

This control constant must be specified for the implicit routines 
and is not used by the other routines. If not specified, the "run” terminates 
with an error message printout. DTIMEI represents a specified time-step and 
is arbitrary, but the governing criterion should be minimum computational 
time with satisfactory tCTiperature accuracy. This means that DTIMEI should 
be specified in conjunction with control constant M.00P ^diich represents the 
maximum number of computational Iterations allowed during each time-step. Since 
each iterative calculation is essentially equivalent to a time-step calcula- 
tion, DTIMEI should be normally greater than NL00P*CSGMIN, where CS(MIN is 
the time-step used in the explicit routines. If savings in computational 
time cannot be met with the same accuracy by using the implicit routines, 
it is more reasonable to use the explicit routines. 

0TIMEL (Minimum Time-Step Allowed) 

This control constant must be specified for subroutine CNFAST and 
is optional for other explicit solution routines. If not specified for 
CHPAST, the "rvm" terminates with an error message printout. DTIMEL repre- 
sents the minimum time-step allowed; for all the ^cpllcit routines except 
CMFAST, if the calculated time-step is less than DTIMEL, the ^*run" terminates 
with an error message printout. For subroutine CNFAST, if the calculated 
time-step of any node, as ^pressed by C^/SG^j and stored in GSGMIN, is less 
than DTIMEL, the tmperature of the nodes not satisfying DTIMEL are calculated 
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using the steady state equations without computational iterations 
(refer to Section 6.3.3 for details on the CNFAST routine). Tlie pur- 
pose of this control constant for CNFASi is to shorten the computational 
time; the danger in its use is that with a large DTIMEL a large number of 
diffusion nodes will receive ths steady state equations without iterations. 
As a restilt, the temperature inaccuracies can be expected to be large. 

DTMPCA (Allowable Diffusion Node Temperature Change) 

This control constant may be optionally specified by the user for 
the Implicit routines and for the explicit routines except CNFAST. DTMPCA 
represents a diffusion— node temperature change criterion between one time- 
step and another. If the maximum diffusion-node temperature change which is 
stored in DTMPCC is greater than DTMPCA, the time-step is shortened to. 

At = .95 * At (DTMPCA/DTMPCC) 

and the diffusion— node and arithmetic— node temperatures re— set to former 
values. The computational procedure is repeated with the smaller time- 
step. DTMPCA serves the same purpose for the diffusion nodes as control 
constant DELXCA provides the arithmetic nodes. 

LAXPAC (Number of Iterations for Linearized Lumped Parameter System) 

LAXFAC is used only In the steady state routine CINDSM and repre- 
sents the number of iterations to be performed on a linear lumped parameter 
system with no updating of elements during a set of LAXFAC iterations. The 
system elements are re-evaluated for the new set of temperatures and in 
turn temperatures are recalculated for another set of LAXFAC iterations 
with a more severe relaxation criterion. The number of iterations will not 
exceed control constant NL00P which represents the total number of itera- 
tions. NL00P will not be met only if relaxation criteria are met during an 
iterative loop and between iterative loops and if the system energy balance 
as stored in BALENG is satisfied (refer to Section 6.5.3 for details). 

NI^0P (Number of Iteration Loops) 

This control constant must be specified for the Implicit and the 
steady state routines; if not specified, the "run" terminates with an 
error message printout. NL00P may be optionally specified for the explicit 
routines since it is used for the arithmetic nodes; if not specified, NL00P 
is set to 1. The value of NL00P to be used depends upon the problem 
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to be solved. For a steady state problem it is not xmusual to have NL00P 
equal to several hundred, whereas for a transient problem the implicit 
routines NL00P should be specified as discussed for control constant 
DTIMEI. In general, a trial and error procedure is required to arrive at 
a suitable value of NL00P. 

OUTPUT (Time Interval for Activating 0UTPUT CALLS) 

. This control constant must be specified for all numerical solution 
routines except steady state routines since the first time-step used is 
generally set to 0UTPUT. The input value is left to the judgment of the 
user. Normally, the output interval is gauged by the length of the run and 
the expected temperature response characteristics. As a "rule-of-thumb" 
the output interval lies between CSGMIN and CSGMAX, with 0UTPUT being 
several times larger than CSGMIN. The values of CSGMIN and CSGMAX can be 
obtained from the output subroutines CSGDMP and RCDUMP.^*** Subroutines 
CSGDMP and RCDUMP are designed to aid in the checkout of thermal problem 
data decks and should be used before making a transient computer run. 

TIMEND (Problem Stop Time) 

The use of this control constant is self-explanatory. For the 
subroutines as they are presently coded, TIMEND must be specified as larger 
than TIME0, otherwise an error message is printed and the "rxm" terminated. 
For the explicit routines j if TIMEND is not larger than TH'IE0 a time-step 
of zero will result and the "TIME STEP TOO SMALL" error message will be 
printed. The implicit routines will print the error message, "TRANSIENT 
TIME NOT SPECIFIED." If a solution is to be terminated by the use of a 
criteria, but the run is not to be terminated, this can be accommodated by 
setting TIMEND=TIME0 when the criteria is met. 

TIME0 ("Old" Time or Problem Start Time) 

This control constant represents the "old" time or the problem 
start time for the transient routines. If not specified, TIME0 is set to 
0.0. An important consideration in the use of TIME0 is that TIME0 may be 
set to negative. 

6.2.4 Time-Step Calculations 

Each numerical solution routine requires the use of a time-step 
that depends upon many considerations, such as the output interval, the end 
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of the problem time, the stability criterion for explicit routines, etc. 

In spite of the xmique solution procedure of each of the numerical solution.' 
routines, the overall time-step calculation procedure for the transient 
routines is essentially identical. The nijmerous time-step checks, as well 
as the selection of the time-step, are indicated below (for definition of 
control constants refer to Section 6.2.3); 

(1) Check that elapsed time, t, does not exceed problem end time. 

If ; TIME0 + 0UTPOT > TIHESD 
Set; 0UTPUT = TIMEND - TIME0 

TIME0 is the bid time 

0UTPUT is the output time interval 

TIMEND is the problem stop time 

(2) Set initial time-step. At, which is stored in DTIMEU (control 

constant for time-step). The initial time step for the SINDA 
numerical routines is as follows; { 

Numerical Routines Initial Time-Step 

EXPLICIT CNFRWD 0UTPUT 

EXPLICIT CNFEDL 0UTPUT 

EXPLICIT CNEXPN 0UTPUT 

EXPLICIT CNDUFE 0UTPUT 

EXPLICIT CNQUIK 0UTPUT 

EXPLICIT CNFAST DTIMEL (minimum time-step allowed) 
IMPLICIT CNMCK DTIMEI (specified time-step) 

DffLICIT CNFWBK DTIMEI 

IMPLICIT CNVAEB DTIMEI 

(3) Check At (stored in DTIMEU) against maximum allowable time-step. 

If ; DTIMEU > DTIMEH 
Set; DTIMEU = DTIMEH 

(4) Check sum of elapsed time since last printout, TSUM, and time- 
step, DTIMEU, against 0UTPUT. 

If; TSUM + DTIMEU > 0UTPUT 
Set; At = 0UTPUT - TSUM 
If; TSUM + At < 0UTPUT 


6 - 23 



and if: TSDM + 2 (At) > 0U1:PDT 

Set : At » 1/2 (OUTPUT - TSUM) 

(5) Store 

Set: DTIMEU = At 

(6) Check DTIMEU against BlnlTimni allowable tinte-step. 

If : DTIMEU < DTIMIL 

Result: An error message is printed and the "run" teanninated 

except for CNFAST, CNBAGK, CNFWBK and CSVAEB. 

(7) Set new time (TIMEN) 

Set: TIMEN = TPRINT + TSDM + At 

TPRINT is the time of the last printout. 

TSUM is the time from the last printout. 

(8) Set mean time (TIMEM) 

Set: TIMEM « 1/2 (TBiS? + TIME0) 

(9) Calculate (or specify) time-step. 

The calculated (or specified) time-step for the SE®A numerical 
routines is as follows: 


Numerical Routines 


Calculated Time-Stet 


EXPLICIT CNFRWD 0.95 * CSGMIN/CSGFAC 

EXPLICIT CNFRDL 0.95 * CSQIIN/CSGFAC 

EXPLICIT CNEOTN 0.95 * CSGMIN * CSGFAC 

EXPLICIT CNDUFR 0.95 * CSOIIN * CSGFAC 

EXPLICIT CNQUIK 0.95 * CSOIIN * CSGFAC 

EXPLICIT CNFAST larger of CSGMIN or DTIMEL 

IMPLICIT CNBACK DTIMEI 

IMPLICIT CNFWBK DTIMEI 

IMPLICIT CNVARB DTIMEI 


CSCMIN “ C^/ZG^^ (minimum value, i = 1,2, . . . ,NND) 
where: C^^ is the capacitance of the ith node 

G^j is the conductance from node i to node j 
CSGFAC is the time-step factor (see above). 
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(10) It should be recognized that individual routines may have slight 
variations to the time-step calculations . 

6.2.5 Computation of Temperatures 

The actual calculation of temperatures, be it for diffusion nodes 
or for arithmetic nodes, represents the end result of a long computational 
procedure with many checks and criteria. Nevertheless, if one confines the 
■discussion to the D0 loops of nodal types, a rather compact but general 
computational pattern becomes apparent. More details are presented in the 
individual sections describing each nximerical solution routine. (Sections 
6.3 - 6.5) 

6. 2. 5.1' Transient Explicit Routines 

For the explicit routines the diffusion and arithmetic nodes are 
treated separately. Diffusion-node temperatures are calculated explicitly, 
whereas the arithmetic-node temperatures are computed implicitly. This 
means that at each time-step an iterative loop is set-up for the arithmetic 
nodes; none is required for the diffusion nodes. 

D iffusion-Node Temperatures 

Calculation of the diffusion-node temperatures follows the 
VARIABLES 1 call^; the computational pattern is: 

D0-L00P (I * 1, NND) On the diffusion nodes is established. 


The functions associated with the variable capacitance C^, the 

variable impressed source q^^, and the variable coefficients (a^^ 
for conduction and radiation), between diffusion-diffusion 

and diffusion-arithmetic nodes are updated at the beginning of each 
time-step. These functional types are described in Section 6. 2. 1.2 
and the computational pattern is indicated in the flow chart of 
Figure 6.2-3. 


Using the updated C^, q^ and the branch heat flow s\im, Qg^» 

and conductance stnn X^, are calculated (refer for example to flow 
chart of Figure 6.3-1). 
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* Variable capacitance (c.) 
impressed source (q^), or 
variable coefficient (G^^) 

Evaluation of Nonlinear Capacitance, 

Source or Conductance 


Figure 6.2-3. 
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Stability criterion C./ E G. . is computed and the smallest value 

is stored in control constant CSGMIN. If CSGMIN 0.0, an error 
message is printed and the "run" terminated. 


Diffusion-node temperatures are calculated by using the appropriate 

finite difference expression associated with each 

routine. These routines and algorithms are Identified as; 

CNFRWD and CNFBDL (Section 6.3.1), uses standard forward- 
difference algorithm. 

CNFAST (Section 6.3.2), uses a modified CNFRI® computational 
procedure to decrease the computational time. 

CNEXPN (Section 6.3.3), uses the exponential prediction method. 

CNDDFR (section 6.3.4), uses DuFort-Frankel method. 

CNQDIK (Section 6.3.5), uses half DuFort-Frankel and half 
exponential prediction metnod. 


Symbolically, the expression for the diffusion-node temperatures may 
be written as. 


T 

i,n+l 




(6.2-3) 


Except for CNFAST the maximum diffusion-node temperature change 
which is stored in DTMPCC is checked against the allowable diffusion node 
temperature change which may be specified by the user via the control 
constant DTMPCA (if not specified DTMPCA = l.OE+8). If DTMPCA is not 
satisfied, the time-step is decreased to. 

At » .95 * At (DTMPCA/DTMPCC) 
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and all temperatures re-set to former values. The computational procedure 
is repeated with the smaller time-step. CNFAST does not allow for the 
recalculation of diffusion-node temperatures, 

Arithmetic-Kode Temperatures 

Calculation of the arithmetic-node temperatures always follows the 
computation of the diffusion-node temperatures and uses "successive point" 
iteration. The computational pattern is as follows; 

Arithmetic-node damping factors DN and DD are established. 

DN = DAMPA (optionally specified user constant, if not specified 
DAMPA = 1.0; factor for the current time-step temperature 
change) 

DD - 1.0 - DN (factor that allows a certain fraction of the "old" 
temperature to be included as part of the temperature change 
for the current time-step) j 

Iterative D0-L00P (K=1,NL00P) is established (NL00P is the number of 

iterations specified by the user, if not specified, NL00P = 1). 


D0-L00P (I*NND, NND + NNA) for the arithmetic nodes is established. 


Impressed source and coefficient (a^^ for conduction and 
aby for radiation) are updated once for each time-step. 


Using the updated and q^, the branch heat flow sum and 
the conductance sum are calculated (refer to flow chart of 
Figure 6.3-2). 
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Arithmetic node temperatures are calculated for each iterative loop 
by using the following "successive point" expression, which is employed 
in all of the routines. 
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( 6 . 2 - 6 ) 
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temperature at kth iteration 

^l.n ^±l.n ^l,k> 

(A « k if j > i and 2 = kt-1 if j < i) 

(a.. and b , . mean updating at time-step, n) 
i3»ii i3»ia 

, b^j = optionally specified (refer to Tables 6.2-1 - 6.2-4) 


DN =: DAMFA (arithmetic node damping factor) 

DD = 1.0 - DN 

The maximum arithmetic-node relaxation temperature change is 
calculated and checked against the allowable arithmetic-node relaxation 
temperature change which may be specified via the control constant ASLXCA. 

This relaxation convergence check is made during each iterative step calcu- 
lation and is used in conjunction with control constant NL00P. Satisfaction 
of either AKLXCA or NL00F during any Iterative step terminates the arithmetic- 
node temperature calculation. 

For each time step, except for CNFAST, the maximum arithmetic-node 
temperature change which is stored in control constant ATMPCC is checked 
against the allowable arithmetic-node temperature change which may be 
specified via the control constant ATMPCA (if not specified, ATMPGA = l.OE+8). 
If ATMPCA is not satisfied, the time-step is decreased to , 

At * . 95 *At (ATMPCA/ ATMPCC) 

and all temperatures re-set to former values. The computational procedure 
is repeated with the smaller time-step. CNFAST does not allow for recalcula- 
tion of arithmetic-node temperatures. 

6. 2. 5. 2 Transient Implicit Routines 

Both diffusion-node and arithmetic-node temperatures are calculated 
by "successive point" iteration. Although these calculations are performed 
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cn the same iterative pass, diffusion node temperatures are evaluated on its 
own computational loop using a specified algorithm associated with a partic- 
ular implicit routine. Calculation of the arithmetic-node temperatures is 
also done on its own computational loop and is identical in all the implicit 
routines. As a matter of fact, arithmetic-node temperatures are calculated 
in the same manner in all the SINDA numerical solution routines. Use of a 
separate computational loop for the diffusion nodes permits the extrapolation 
of diffusion-node temperatures provided acceleration of convergence criterion 
is met (refer to Section 6.2.7). 


Diffusion-Node Temperatures 

In order- to facilitate the discussion to follow on the computa- 
tional procedure, it is convenient to examine the forward-backward finite 
difference expression. 


13 


X At 


e T- - + (1 - 3) T, 1 , (6.2-7) 

forward backward 


where; 3 = factor with range 0 < 3 £ 1/2 

P P 4 4 

T. j = q. + E a.. (T. - T. ) + S Ob.. (T7 -T^ ) (6.2-8) 

forward ^x,n lj ,n a ,n x,n^ i3 ,n 3 


\ackward ^i,n ®ij ,n^^j ,fc+l"^i,k+l^'‘’^f^ ^^ij ,kfl“^i,k+l^ 


1 s 1,2, . . »,N 

^j,n* ’^j.kfl “ <^°“®tant, N < i < p 
n ® nth time-step ; k = kth iteration within a given time-step. 
^i*^i*^ij*^ij ” optionally specified (refer to Tables 6.2-1 — 6.2-4) 

Any value of 3 less than one yields an implicit set of equations 
which must be solved simultaneously. For values of 3 less than or equal to 
one-half equation (6.2-7) represents an unconditionally stable set of equa- 
tions, whereas values of 3 greater than one-half yields a set of equations 
with conditional stability. 

The standard implicit algorithm used in subroutine CNBACK follows 
directly from equation (6.2—7) by letting 3 = 0, whereas the Crank-Nicolson 
method used in subroutine CNFWBK follows by letting 3 “ 1/2. Subroutine 
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CNVARB uses a variable factor which is based upon the ratio of CSQIIN/DTIMEUj 
this ratio is internally calculated in CNVAEB (refer to Section 6. 4. 3. 2). 

In order to simplify the presentation, the following notation is used. 

lor CNBACK (3 = 0): 


Qi 


= q. + C. T. 
x,n i,n x,n 


( 6 . 2 - 10 ) 
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from ith node, called radiation damping (refer to 
Section 6.2.6 for details) 

= 0,if radiation is not present 

For CNFWBK (3 — (note equation (6.2-7) is multiplied by 2); 
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= same as equation (6.2-11) 
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(6.2-15) 


(6.2-16) 


G . . » same as equation (6.2-13) 

^^i^ave ~ equation (6.2-14) 

For CNVARB (variable 3') (note that equation (6.2-7) is multiplied by 2, 
so that 3* = 23 now ranges, 0 ^ B’ 1.0): 
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loss from 1 th node, called radiation damping (refer to 
Section 6.2.6 for details) 


0 , if radiation is not present 


i » 1,2,..., N 

3 = 2.0*CSGMIN/DTIMEU (range allowed, 0 < 3* < 1.0) 

T. i T. , - constant, N < j < p ' (p is the total number of nodes) 
3,a 3»*^ — 


( 6 . 2 - 20 ) 


n = nth time-step; k = kth iteration 


Ci* q^, a^^, b^j 




/At 


may be optionally specified (refer to Tables 6.2-1 - 6.2-4) 


Calculation of the diffusion-node temperatures follows VARIABLES 1 
call; the computational pattern is: 

Iterative D0-L00P (kl=l,NL00P) for the total nodal system is established. 
First Iterative Loop: 

DO-LOOP (I=1,1]M}) on diffusion nodes is established. 

The functions associated with the variable capacitance C^, the 
variable impressed source q^, and the variable coefficients Gj^ (a^^ 
for conduction and for radiation) between diffusion-diffusion 

and diffusion-arithmetic nodes are updated once for each time-step. 

These functional types are described in Section 6 . 2. 1.2 and the com- 
putational pattern is indicated in the flow chart of Figure 6.2-3. 
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Ail known quantities (those evaltiated at time-step n) are summed and 
are identified by the symbol (equations 6.2-10, 6.2-15 and 6.2-17). 

CSGMIN is evaluated. 

Radiation damping is used; average radiation heat loss, from 

the ith node is evalxiated (refer to Section 6.2.6). 

For CNVARB, B’= 2.0*CSGMIN/DTIMEU is calculated. 

The diffusion-node temperatures are calculated by ’’successive point” 
iteration (actually CNBACK and GNFWBK have slightly different first 
iterative pattern than CNVARB but the difference is hot significant). 

DN » DAMPD (use specified diffusion node damping factor, 
if not specified, DAMPD = 1.0) 


DD » 1.0 - DN 

1 

For CNVARB, the diffusion-node relaxation temperature change is 
calculated; maximum value is stored in DRLXCC. 

Second and Sueceeding Iterative Loops ; 

With the iterative loops after the first, those quantities C^, q^, and 
6^ which were updated during the first iteration are held constant. 

Diffusion-node temperatures are found by using equation (6.2-21). 


The diffusion-node relaxation temperature change is calculated and the 
maximum value stored in DRLXCC. 


Check of DRLXCC against DRLXCA (allowable maximum diffusion-node relaxa- 
tion temperature change) is made after the arithmetic-node temperature 
calculations. 

Each third iteration, a check on solution convergence is made; if con- 
vergence is occurring linear extrapolation to accelerate convergence 
is made (refer to Section 6.2.7). 

Arithmetic-Node Temperatures (if any) 


During the first iterative loop the impressed source q^ and 
coefficient G, (a,, for conduction and ab.. for radiation) between arithmetic- 

K Xj ij 

arithmetic nodes are updated once each time-step. On every loop, arithmetic- 
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node temperatures are calculated using "successive point" iteration. The 
finite difference algorithm is presented in Section 6, 2. 5.1 (equation 6.2-6). 

The arithmetic-node relaxation temperature change is calculated 
and the maximum is stored in .ABLXCC. 

During Each Iterative Loop After the First 

Both DEI^CC and AKLXCC are checked against DRUCCA and ARLXCA, 
respectively. If both DRLXCA and ARLXCA are satisfied, the iteration 
ceases . 

If L00PCT equals NL00P the message "REIAXATION CRITERIA NOT MET" 
is printed. 

Both the calculated maximum diffusion-node and arithmetic-node 
temperature change (stored in DTMPCC and ATMPCC, respectively) are checked 
against the corresponding allowable temperature change stored in DTMPCA 
and ATMPCA. If DTMPCA is not satisfied, the time-step is decreased to. 

At = .95*At(DTMPCA/DTMPCC) 

and all temperatures re-set to former values. The computational procedure 
is repeated with the smaller time-step. 

If ATMPCA is not satisfied, the time-step is decreased to. 

At = .95*At (ATMPCA/ ATMPCC) 

and all temperatures re-set to former values. The computational procedure 
is repeated with the smaller time-step. 

6. 2. 5. 3 Steady State Routines 

Diffusion nodes and arithmetic nodes are treated separately in 
CINDSS and CINDSL even though from a physical standpoint a distinction 
between diffusion nodes (nodes with capacitance) and arithmetic nodes (nodes 
with no capacitance) doesn't exist. Thus, the set of control constants for 
the diffusion nodes and another set of control constants for arithmetic 
nodes are similar to those used in the transient routines. No distinction 
in the type of nodes is made in CINDSM. 

The computational procedure to be discussed applies only to CINDSS 
and CINDSL:' CINDSM is considerably different (refer to Section 6.5.3). 
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Dlffusion-Kode Tesiperatures (iiodas specified with capacitance even though 
the problem is steady state) 

Aa iterative D0-L00P (K1-1,NL00P) is established. 

Within this iterative loop a D0-L00P (1=1, MD) on the diffusion 

nodes is made. The functions associated with the impressed source and the 

variable coefficients G, (a., for conduction and 0 b.. for radiation) between 

k xj 

diffusion-diffusion and diffusion-arithmetic nodes are updated each iteration. 

Diffusion-node temperatures are calculated using "block” iteration 
for CINDSS and "successive point" iteration for CIHDSL. 


"Block" iteration (CINDSS) i 


T 4 . t . , = DD*T . . + DN* 
i,k+l x,k 


(^i»k ^ ^ij ,k ^j ,k) 
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j=l 


( 6 , 2 - 22 ) 


ij,k 


®ij,k “ “lj,k ^ ^^ij,k ^^j,k t ^i,k^^^j,k ^i,k^ 

DN = DAMPD (diffusion-node damping factor) 

DD = 1.0 - DN 

i = 1,2,...,NND (number of diffusion nodes) 

k * kth Iteration; p = total number of nodes 

‘*i*^ij’^ij “ optionally specified to Tables (6.2-1 - 6.2-4) 

T. - constant, (NND + NNA) < j p (NNA is the number 

of arithmetic nodes 

"Successive point" iteration (CINDSL) ; 






i=i-H 


, k ,k) 


P 
2 G 
j=l 


(6.2-23) 


ij,k 


=«,k " ^j.k * '«'l 3 ,k + \.k> 


(Jl» k if j > i and £ » k +1 if j < i) 


DN = DAMPD 

DD = 1.0 - DN 


6 - 35 



1 = 1, 2, .... (NND + KNA) 

Ic »* kth iteration; p = total nimber of nodes 

“ optionally specified to Tables (6.2-1 - 6.2-4) 
T. , “ constant y (!?HD + iNNA) < j ^ P (NNA is the total 

J j lc 

number of arithmetic nodes) 

Diffusion-node relaxation temperature change is calculated and 
the maximum is stored in DBLXCC. 

Arithmetic-Node Temperatures (nodes specified with no capacitance) 

Within this iterative D0-L00P a D0-L00P (I=NND+1, NND + NNA) is 
established. 

The functions associated with impressed source and variable 
coefficients ^^ij conduction and b^^ for radiation) between 
arithmetic- arithmetic nodes are updated each iteration. 

Arithmetic-node temperatures are calculated using "successive 
point" iteration. 




fl.k 

+ AN* — 


P P 

2 G . . 1 ^ T . ■t,-x "b 2 6 . . 

1=1^*^ 3,k-fl 


,k ^j,k) 


j 


\ ®ij.k 


®ij,k * ®lj,k'‘'^\j,k ^ ^i,P 

(A * k if j _> 1 and A = fcfl if j < i) 

AN = DAMPA (arithmetic-node damping factor) 

AD » 1.0 - AN 

i = (NNDfl) , (NND+2 ) , ...» (NND + NNA) (number of 

arithmetic nodes) 

k kth iteration 

p total number of nodes 

T . . = constant, (NND + NNA) < j p 

J f 


The arithmetic-node relaxation temperature change is calculated 
and the maximum value is stored in ARLXCC. 


.2-24) 
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During Each Iterative Loop 


Both DRLXCC and ASLXCC are checked against DRLXCA and ABLXCA, 
respectively. If both relaxation criteria, DRLXCA and ARLXCA, are satisfied, 
the iteration ceases. 

If both relaxation criteria, DRLXCA and ARLXCA, are not met with 
NL00P iterations, the message "ITERATION COUNT EXCEEDED, NL00P * " is 

printed. 

Energy balance of the system is calculated and is stored in con- 
trol constant ENGBAL. 


6.2.6 Radiation Damping 

Radiation damping denotes an averaging of radiation heat loss 
technique used to prevent or minimize large temperature oscillations . This 
method is currently employed in only the implicit routines. The technique which 
is original with J. D. Gaskl is based upon practical and computational con- 
siderations, Solution of numerous problems without large temperature 
oscillations indicates the effectiveness of the approach. 


The radiation averaging technique is relatively simple conceptually 
and rather easily incorporated in the numerical solution routines . The com- 
putational pattern is such that the diffusion nodes are encountered sequen- 
tially. Let the encountered node be the 1th node. A check is made for the 
presence of a radiation Goefflcient, Gj^ ith node. If one or 

more radiation connections is present, the radiation heat loss , 


the 1th node is calculated based upon the previous temperature T 


i,k* 


^‘*i^rl " f '^^ij,n ’^i,k 

J* 


(6.2-24) 


where, j = all radiation connections to node i 
n * nth time-step 
k * kth iteration 


Oslng a second temperature (T^ y}2* found as follows : 

"l.k>2 - 
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where. 


Q =C. T, "f" H a.. T, , , - + E a . . T , 

^sum X i,n ^x,n x^.n j ij »n j »k 


+ Z 0b.. + Z 0b., - 


(6.2-26) 


G = C, + Z a., 
sum X i3,n 


(6.2-27) 


Note that in the evaluation of (T. .).,, the damping factor DAMPD is not used. 

Note further that does not contain since it is accotmted for in 

the radiation loss term, ^ 

Now a second radiation heat loss based on (T. . )_ is found, 

X y <6 


^‘*i^r2 ~ j ^^ij,n ^’’^i,k^2 


(6.2-28) 


Equations (6.2-24) and (6.2-28) are then averaged, 
^%^ave “ ^^‘^i^rl ^ 


(6.2-29) 


This average radiation heat loss from an ith node is used in the 
diffusion-node finite difference algorithm as follows. 




^^sxjm ^*^i^ave^ 


(6.2-30) 


where, DN = DAMPD 

DD = 1.0 - DN 

^^i^ave " average radiation heat loss (equation 6.2-29) 

G “ C . + Z a . . 

sum i ^ ij ,n 

^svun " form shown by equation (6.2-26). The actual 

expression depends upon algorithm. Equation (6.2-26) 
is for the standard implicit method. 


The reason behind the use of (q.) is that if the initial 

^x ave 

temperature T^ is too large, the heat loss from the ith node, ( 9 £)j -2 
would then be too large. As a result the evaluation of (T^ with 
(qi)ri would yield a temperature that is too low. Thus, the averaging of 
of Slid ^'^i^r2 much closer to the true heat loss from the 
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ith node. If ^ is too sraall then (T^ would be too large; the 
averaging scheme still holds true. 

6.2.7 Acceleration of Convergence by Extrapolation Technique 

Several of the SI1®A ntimerical solution routines use an extrapola- 
tion technique to accelerate convergence of the iterative procedure. The 
extrapolation technique is used in the implicit routines CNBACK, GNFWBK, 
and CNVARB for the iterative tenperature solution of the diffusion nodes, 
but is not used for the iterative temperature solutions of the arithmetic 
nodes. The extrapolation method is also used in the steady state routines 
CINDSL and CINDSM for the iterative temperature solution of both the 
diffusion and the arithmetic nodes. 


6. 2. 7.1 Extrapolation Technique 

The extrapolation is baaed on a zero temperature difference con- 
dition which is defined to be a point where the temperature change of a 
particular node over two successive Iterations is zero. The governing 
equations are developed as follows: 

Consider the temperatures of an ith node at three successive 
iterations as shown in Figure 6.2-4a. Let these temperatures, which are 
as sinned to be successively decreasing (or increasing), be denoted as. 


^i,k-2* ’^i,k-l ^“‘*/*^i,k 


where, k is the present iteration 

k-1 is the previous iteration 

k-2 is two iterations before the kth iteration 


By taking the differences. 


^^i,k-l \,k-2 “ ^i,k-l 




and plotting these temperature differences as a function of iterations, 
the iterative point of zero temperature difference can be found by linear 
extrapolation as shown in Figure 6. 2-4b. The corresponding expression for 
the line is found by using the point, AT^ at I = k and the slope, 

(AT^^j^ - / (k-(k-D), to yield,’ 
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k-2 k-1 k 

No. of Iterations, I 

Figure 6.2-4a. Temperature (ith) vs. No. of Iterations 



k-1 k R 

No. of Iterations, I 

Figure 6.2-4b. Temperature Difference vs- No. of Iterations 



k-X k K 

No. of Iterations, I 

Figure 6.2-4c. Extrapolation of Temperature (ith) to New Value 
Figure 6.2-4. Method of Extrapolation to Accelerate Convergence 
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(6.2-31) 


where. 


* “i.k + <"i,k - « - « 

I =» iterations 


Since at the zero temperature difference condition, AT. _ “0, the expression 

X j X 

for the extrapolated iterations, >= (K - k), is found to he, 

^*-“i.k/<“i.k-“i,k-i> 

Now, by extrapolating the line established by the temperatures, T^ and 
y to the line I = K, as shown in Figure 6. 2-4c, the extrapolated tempera- 
ture T. ^ is found. Tlie expression is readily found to be, 

^1,1 - ^i.k + «l.k - \,k-l>« - « 

Since I ’= K and K - k = K^, equation (6.2-33) becomes. 


^l,K ■ ’^i.k '1 + V - \ ^l,k-l 


(6.2-34) 


6.2, 7.2 Frogramming Considerations 

Each applicable node is tested at the completion of each third 
iteration to determine if the extrapolation method should be applied. If 
is calculated to be less than or equal to zero, extrapolation is neglected 
since the error function is diverging. If is calculated to be greater 
than zero, a new temperature is calculated based on equation (6.2-34) j how- 
ever, to avoid problems associated with a nearly-zero slope of the line 
representing the temperature difference vs, number of iterations relation- 
ship (Figure 6.2-4b), K^, is set to a number K^; otherwise, could be a 
very large ntimber. For the implicit routines, CNBACK, CNFWBK, and CNVARB, 

K “ 10. For the steady state routine CINDSL K = 8 and for steady state 

IB IB 

routine CINDSM a criterion based upon the maximum temperature is used. 


6. 2.7. 3 Routines Using Acceleration of Convergence 

SINDA numerical solution routines that employ the acceleration 
of convergence features are: 

CINDSL, CINDSM Steady state routines 

CNBACK, CNFWBK, CNVARB Tr^slent implicit routines 

6.2. 7.4 Comment on Acceleration of Convergence 

Neither an extensive study on the value of the acceleration 
convergence feature has been made, nor has one been reported, but the 
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, limited results presently available indicate that for the steady state 
routine CIM>SL, the number of iterations is reduced approximately 20%. 
Results are not available for the implicit routines. 

A study of the acceleration of convergence feature is made 
difficult because the method is not a user option in the applicable SIMDA 
numerical solution routines. Thus, the user must be sufficiently versed 
with the routines in order to delete the acceleration of convergence 
feature. 

6.2.8 Other Characteristics of the SISDA Numerical Solution Routines 

6. 2. 8.1 Units 

SIRDA, as presently coded, requires that the temperatures must be 
specified in degrees Fahrenheit (®F) since the conversion factor to obtain 
degrees absolute is internally set at 460.0. This means that the units 
must be consistent with “F (or ®R). The execution routines as presently , 

I 

• coded do not permit the use of other units. . ' 

6. 2. 8. 2 General Comments on Computational Features 

Many of the computational features such as radiation damping are 
original with J. D. Gaski. No theoretical proofs are offered since a 
practical "gut-feel" development was often used in lieu of a sophisticated 
mathematical approach; the features, in general, appear to meet the 
intended objectives. It should be particularly noted that the numerical 
solution routines are computationally similar; within a particular numerical 
solution class explicit, implicit or steady state, the computational 
similarity is even more pronounced. Yet on the other hand, similarity of 
patterns are brokaa for no particular reason other than the programmer's 
whim. 
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6.3 


Routines 



SINDA explicit solution routine number six. These are identified 
as follows : 

CNFRWD Conditionally stable explicit forward difference. 

Requires short pseudo-compute sequence (SPCS) . 

CNFRDL Identical to CNFRWD except that the long pseudo- 
compute sequence (liPCS) is required. 

CNFAST Modified CNFRWD for accelerated forward differencing. 
Requires short pseudo-compute sequence (SPCS) . 

- CNEXPN Unconditionally stable explicit differencing using 
exponential prediction. 

Requires short pseudo-compute sequence (SPCS). 

CNDUFR Unconditionally stable explicit differencing using 

I 

DuFort-Frankel method. 

Requires short pseudo-compute sequence (SPCS). 

CNQUIK Unconditionally stable explicit differencing using 
a combination of half CNEXPN and half CNDUFR. 

Requires short pseudo-compute sequence (SPCS). 

A detailed description of each explicit routine is presented on 
the pages to follow with heavy reliance upon the general description of 
Section 6.2. A brief description of these routines is summarized first. 

CNFRWD uses an explicit forward differencing algorithm and 
requires the short pseudo-compute sequence (SPCS). The explicit method 
is characterized by computational simplicity and stability limitations. 
Since the allowable time-step is governed by the smallest time constant 
of the network, care must be given in reducing the physical system to a 
reasonable lumped-parameter model. Arithmetic-node temperatures are 
calculated by "successive point" iteration. 

CNFRDL is identical to CNFRWD except that CNFRDL requires the 
long pseudo-compute sequence instead of the short pseudo compute sequence. 
CNFRDL requires slightly less solution time than CNFRWD but the difference 
is not significant; CNFRDL does require more core storage, however. 
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CNFAST represents a modified CN¥RVJD with the modifications 
intended to decrease the computational time. A user specified control 
constant DTIMEL which contains the minimum time-step allowed is used as a 
criterion for isolating those diffusion nodes that are to receive the 
steady state calculations. A large pocket of internally converted diffusion 
nodes can present considerable accuracy problems. 

CNEXPH uses an unconditionally stable explicit method with the 
Intent to reduce computational time at the expense of temperature accuracy. 
If accuracy is an important consideration, another routine such as CNTRWD 
would be a better choice. As a note of Interest, CNEXPN solutions tend to 
lag in time the true solutions, 

CNDUFR uses the unconditionally stable DuFort-Frankel method 
with the intent to reduce computational time by using time-steps greater 
than those allowed with the conditionally stable explicit methods. Again 
accuracy may be compromised. CNDUFR solutions tend to lead in time the i 
time solutions, 

CNQUIK uses half CNEXPN and half CNDUFR. Why? Since CNEXPN 
solutions tend to lag in time and CNDUFR solutions tead to lead in time, a 
combination may yield better solutions. Preliminary results indicate that 
CNQUIK solutions are more accurate than either CNEXPN or CNDUFR for the 
same computational time. 
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6.3.1 Subroutines; CNFRWD and CHFBDL 

6. 3. 1.1 General Comments 

Subroutines CNTBWD and CNFBDL are numerical solution routines 
that use the forward finite difference explicit approximation^^* of the 
parabolic differential equation. CNFRWD and CNFEDL are identical except 
that CNFRDL requires the short pseudo compute sequence (SBCS) whereas 
•CNFRDL requires the long pseudo compute sequence (LPCS) . The need for 
both routines becomes apparent when it is understood that if a steady 
state numerical solution routine is followed by a transient numerical 
solution routine, both routines must have consistent PCS (LPCS or SPCS) . 

As a note of interest, each arithmetic node receives the long pseudo 
compute sequence (LPCS) but this is done internally by the program. 

The forward finite difference explicit method as used in CNFRWD 
and CNFRDL is the conventional Euler method that neither provides a check 
on the accuracy nor does it provide any scheme of correction once the 
temperature values are calculated except for the arithmetic nodes which 
are reiterated NL00P-times, The explicit method is characterized by com- 
putational simplicity and stability limitations with the temperature 
error at any time point being of the order At, 0(At), provided the stability 
criterion is satisfied. For a rapidly .changing boundary condition, such 
as a heat source^ there is no assurance that the calculated temperatures 
are, accurate during the transient period, particularly near the start of 
the transient, even though the stability criterion is satisfied. Since 
the allowable time step is governed by the smallest time constant of the 
network, care must be given in reducing the physical system to a lumped- 
parameter model. Nonlinearity due to the presence of thermal radiation 
exchange or temperature- time varying coefficients can lead to numerical 
solution difficulties j the presence of arithmetic nodes can also present 
difficulties. These routines offer a number of control constants many 
of which can be optionally specified by the user to affect the numerical 
results . 

Even with the experience gained through the use of these 
routines, no realistic criteria can be stated except for the qualitative 
guidelines indicated above. It is thus recommended that the user becomes 
familiar with various control constants and their role. The presentation 


6 - 45 



to follow is inteaded to provide the instructional information. 

6.3. 1.2 Finite Difference Approximation and Computational Algorithm 

The forward finite difference explicit formulation of the 
lumped parameter heat balance equations was presented in Section 5.2.1. 
For convenience) the expression is repeated here. 


(T . - T ) 
^i At 


H,n 





i,n 


P 

+ . S 
3=1 


Ob 


ij 




(From equation 5.2-1 of Section 5.2.1) 


where, i = 1,2,...,N 

T = constant, N < j < p 
j ,n 

p *= total number of nodes 
At “ time-step 
n = nth time-step 


By letting ^ 

(5.2-1) becomes. 


a., „ + Ob.. „ (T^ „ + )(T. ^ 

ij ,n ij ,n 3*^ i»n 3*®- 


^) , eqtiation 


(T. - T . ) 

i,n-M x.n 
At 


P 

= q. + 2 G.. 
3 -»n X3 ,n 


3 i,n 


(6.3-1) 


The algorithm as used in the subroutines for the diffusion nodes and for 
the arithmetic nodes may be expressed as follows. 


Diffusion Nodes 


= T, 


‘i,n+l i,n 


C. 


[ *^i,n 


P 

2 

j=l 


’ij.n 


■ i J 

i,n 


] 


(6.3-2) 


where, n * nth time-step 

At = time-step (refer to Section 6,2.4) 

i * 1,2,...,NND (number of diffusion nodes) 

T. * constant, (NND + NM) < J < p (NNA is the number of 
J»n . 

arithmetic nodes and p is the total number of nodes) 
C^, q^, a^j , b^j ® may be optionally specified (refer to 

Tables 6.2-1 through 6.2-4). 
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Arithmetic Nodes (if any) 


^i,k+l “ ^i,k 


(q. + Z G.. T. 

\ jsi 13 ,n 3» 


1(^1 + Z G^. 
i*l+l ^ 


,n ,k) 


P 

Z G 
4“l 


(6.3-3) 


ij,n 


where, k = kth iteration loop; i * (NND + 1), (Kl® 4- 2),.. .,(NND + NNA) 

C^,q^,a^^,b^j = optionally specified (refer to Tables 6.2-1 - 6,2-4) 

T. , * constant, (NHD + MA) < j < p (NNA is the number of 
3 » K 

arithmetic nodes and p is the total number of nodes) 

DN = DAMPA (arithmetic node damping factor, refer to Section 6. 2. 3. 2) 
- DD » 1.0 - DN 

G.. = a.. + ab.. (T? . + T? ,)(T. „ + T. ,) 

i3»n i3»n ij ,n '3,2- i,k'' x,k' 


(A * k, if j > i and 2 = k+1, if j < i) 


6. 3. 1.3 Comments on the Computational Procedure 

The important steps of the computational procedure used in 
subroutines CNFRWD and CNPEDL are indicated in Table 6.3-1. For a 
detailed step-by-step computational description, the user must examine the 
computer listings for CNFRWD and CNFEDi presented in Appendix A, but some 
general computational details are given in Section 6. 2. 5.1. Both CNFRWD 
and CNFRDL use essentially the same computational steps with the difference 
occurring in the calculation of the diffusion-node temperatures as shown 
in the flow chart of Figure 6,3-1; a flow chart for the calculation of the 
arithmetic-node temperatures is shown in Figure 6.3-2, A functional flow 
chart of CNFRWD and CNFRDL is shown in Figure 6.3-3. The difference 
between CNFRWD and CNFRDL is due to the use of the short pseudo-compute 
sequence (SPCS) by CNFRW) and the use of the long pseudo-compute sequence 
(LPCS) by CNFRDL. 

All diffusion-node temperatures are calculated by a two-pass 
operation prior to the calculation of the arithmetic node temperatures. 

On the first pass the psevido-compute sequence for the diffusion nodes is 
addressed and the heat flow is calculated and the direction determined for 
each conductor encountered; the appropriate heat flow and conductance 
summations are performed. Refer to Section 6. 2.5.1 for more details on 
the computational procedure. 
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The stability criterion of each diffusion node is calculated 

"I and the minionim value is placed in control constant CSQ!IN. The time-step 
used (stored in control constant DXIMEU) is calculated as 95% of CSGMIN 
divided by control constant CSGFAC which is set at 1.0 unless specified 
larger by the user. A "look ahead" feature is used when DTIMHI is calculated. 
If one time-step will pass the output time point the time-step is set to 
lie on the output time point; if two time-steps will pass the output time 
point, the. time-step is set so that the end of the two time-steps will lie 
on the output time point. DTIMEO is checked against both DTIMEH and DTIMEL. 

If DTIMEO exceeds DTIMEH, DTIMEU is set equal to DTIMEH, and if DTIMEU is 
less than DTIMEL, the "run" is terminated. DTIMEL is internally set to 
zero if not specified and DTIMEH is set to l.OE+8 if not specified. The 
maximum diffusion node temperature change over a time-step is placed in 
control constant DTMPCC and is checked against the allowable diffusion node 
temperature change stored in the optionally user specified control constant 
DTMPCA which is not Specified is set to l.OE+8. If DTMPCC is larger than 
DTMPCA, DTIMEU is shortened and the calculations repeated. Refer to 
Section 6.2.4 for detailed procedure on time-step calculation. 

) The user may iterate the arithmetic node calculations during a 

time-step by specifying control constant NL00P and adjust the solution by 
the use of ARLXCA. The maximum arithmetic node temperature change over an 
iteration is placed in control constant ARLXCC and is checked against the 
arithmetic node temperature change criterion stored in ARLXCA. Satisfaction 
of either NL00P or ARLXCA terminates the iterative process for that time-step. 
If the arithmetic node iteration count exceeds ML00P the results are retained 
and computation proceeds without user notification . The maximum arithmetic 
node temperature change over the time^step is stored in control constant 
ATMPCC and is checked against the allowable temperature change stored in 
ATMPCA. If larger, the time-step is shortened and the calculation repeated. 
The user may also specify the control constant DAMPA in order to dampen 
possible oscillation due to nonlinearities. 

6. 3. 1.4 Control Constants 

Control constants 0UTPUT and TIMEND (> TIME0) must be specified 
as indicated in Table 6.2-5 and described in Section 6. 2. 3. 2; otherwise the 
i "run" will terminate with an error message. The function of optionally 
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specified control constants AELXCA, ATMPCA, BACKUP, CSGFAC, DAMPA, DTIMEH, 
DTIMEL, DTMPCA, NL00P, and TIME0 is described in Section 6. 2.3. 2. Note 
particularly that TIME0 may be set negative and that NL00P is set to one 
if not specified. 

6. 3. 1.5 Error and Other Messages 

If control constants 0DTPUT and TIMEND are not specified, the 
following error message will be printed for each, 

0UTPUT "N0 0UTPUT INTERVAL" 

TIMEND "TIME STEP T00 SMALL" 

The reason for the TIMEND error message is that a direct check on TIMEND 
is not made; the resultant error message just happens to be a quirk in 
the coding. 

If the short pseudo-compute sequence SPCS is not specified, 
the error message will be, 

"CNFRWD REQUIRES SHORT PSEUD0-COMPUTE SEQUENCE" 

If the long pseudo-compute LPCS is not specified, the error 
message will be, 

"CNFRDL REQUIRES LONG PSEUDO-COMPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficient 
(NDIM < (NND + NNA) ) , the message will be , 

" L0CATIONS AVAILABLE" 

Note that the number printed will be negative Indicating the additional 
storage locations required. 

If the time-step used is less than the time-step allowed 
(DTIMEL) which may be optionally specified by the user, the message will be, 

"TIME STEP T00 SMALL" 

If GSGMIN < 0, the message printed will be, 

"CSGMIN ZER0 or NEGATIVE" 

Checks on the control constants , the pseudo-compute sequence and 
the dynamic storage allocation are made in the following sequence with the 
"run** terminating if a single check is not satisfied. 
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0UTPUT, pseudo-compute sequence, dynamic storage locations 

It should be particu3.arly noted that no message is printed if 
AELXCA is not satisfied with NL00P iterations; ARIICCA and KL00P are 
optionally specified control constants. 
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1 . 


2 , 

3. 

4. 

5. 

6 . 
7. 


8 . 

9. 

10 . 

11 . 


12 . 

13. 

14. 


Table 6.3-1. Basic Computational Steps for CNFRWD and CNFRDL 

Specification of control constants (all control constants are pre-set to 
zero). Control constants 0UTPUT and TIMEND must be specified. SPCS is 
required for CNFE.WD and LPCS for CNFRDL. (Refer to Table 6.2-4 for 
nominal values and Section 6.2. 3. 2 for description.) 

Sufficiency check on dynamic storage. Requirements NND + NNA (NND = 
diffusion nodes and NNA - arithmetic nodes) , 

Setting and/or calctilation of time-step* At. (Refer to Section 6.2.4 
for detailed procedure. ) 


Setting of source and diffusion node dynamic storage locations at zero. 
Calling of VARIABLES 1. (Refer to Section 6. 2, 2. 2.) 

Cheeking of BACKUP. (Befer to Section 6. 2. 3. 2.) 

Calculation of diffusion-node temperatures. (Refer to Section 6. 2. 5.1 
for description and to flow chart of Figure 6.3-1.) 

Diffusion-node temperatures are calculated by using: (refer to Section 
6.3.1. 2.) 


i,n+l 


T, + 
i,n 


AT 


i,n 


where. 




[ 


^i,n 


P 

+ Z 
j»l 


G.. (T. 

i3,n ' 3 ,n 


r.5 

i,n_ 


Erasure of all temperature calculations for latest time-step if allow- 
able temperature change criterion DTMPCA is not satisfied and recalcula- 
tion of temperatures with reduced time-step. 


Calculation of arithmetic-node temperatures; if the number of iterations 
equals NL00P the temperatures are retained without user modification. 

(Refer to Section 6. 2. 5.1 for description and to flow chart of Figure 6.3.2) 


Erasure of all temperature calculations for latest time-step if allowable 
temperature change criterion ATMPCA is not satisfied and recalculation 
of temperatures with reduced time-step. 

Setting of BACKUP to 0.0 and the calling of VARIABLES 2. 

If BACKUP is nonzero, temperatures are re-set to former values and the 
computational procedure repeated. 

Advancing of time, cheeking of time to print, and the printing at the 
output interval. 

Calling of 0UTPUT CALLS. 

Checking for problem end-time stored in user specified control constant 
TIMEND. 
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Figure 6.3—1. QSUM and GSUM for "Block" Diffusion— Node 
Temperature Calculation, CNFRWD and CNFRDL 
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Figure 6.3-3. . Functional Flow Chart for CNFRWD and CNFRDL 
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6.3.2 


Subroutine; CNFAST 


6. 3. 2.1 General Coimnents 

Subroutine CNFAST, which requires the short pseudo compute sequence 
(SPCS) represents a modified CNFRWD with the modifications intended to 
decrease the computational time. Use of CNFAST requires a user specifica- 
tion of control constant DTIMEL which represents the minimum time-step 
allowed in addition to control constant 0UTPUT, With minimum computational 
time and adequate temperature values as the objective, the computational 
procedure is simplified. A number of checks on control constants are 
eliminated and temperature nodes with CSGMIN less than the allowable time- 
step, DTIMEL, are calculated using the steady state equations. 

Although experience on the use of CNPAST is rather limited at this 
time, it is clear that the user specified DTIMEL should be sufficiently 
small that only a small number of the diffusion nodes should receive the 
steady state equations. These steady state equations are computed only once 
during a time-step and thus are not treated computationally the same as che 
other user-specified arithmetic nodes. A large pocket of internally con- 
verted diffusion nodes would lead to large temperature inaccuracies. 

6. 3. 2. 2 Finite Difference Approximation and Computational Algorithm 

The finite difference expressions for CNFAST are the same as those 
indicated in Section 6. 3.1.2 for subroutines CNFRWD and CNFRDL, but the 
application of these equations in the computation procedure is different. 

Diffusion Nodes 


If the user specified control constant, DTIMEL, which represents 
the maximum time-step allowed as specified by the user is less than or equal 
to CSGMIN, the diffusion node temperature is calculated as. 


J.1 “ T. + 

i,n+l i,n 
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— q. 
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■ij,n ^^j,n 
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6.3-4) 


where. At 

^i*V®ij »^ij 
i 


= time-step (refer to Section 6.2.4); n = nth time-step 
® optionally specified (refer to Tables 6.2-1 - 6.2-4) 

= 1,2,...,NND (of diffusion nodes with DTIMEL < CSGMIN) 


j.n 


'ij.n 


constant, (NND + NNA) < j ^ p (NNA is the number of 

arithmetic nodes and p is the total nvimber of nodes) 

a + Ob.. ^ (T^ ^ + T? „)(T. „ + T. j 
rj ,n ij j i»^^ J i*'^ 
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If DTIMEL > CSQilN, the time-step is set at DTIMEL and the 
diffusion node temperature calculated with no iterations as. 


T 

■^i,n+l 






P 

3“1 





(6.3-5) 


where, n means the nth time-step 

i * 1,2,...,NM) (number of diffusion nodes with DTIHBL > CSGMIN) 


Arithmetic Nodes (if any) 

The arithmetic-node temperatures are calculated in the same manner 
as in CNFRWD (Section 6. 3. 1.2) or refer to Section 5.2.3 for finite difference 
algorithm. 


6. 3. 2. 3 Comments on the Computational Procedure 


The Important steps of the computational procedure used in 
subroutine C3SJFAST are indicated in Table 6.3—2 and a functional flow chart 
is shown in Figure 6.3-4. For a detailed computational description, the 
user should examine the computer listing for GNFAST in Appendix A, but some 
general computational details are presented in Section 6. 2. 5.1. The compu- 
tational procedure is similar to the one used in CNFRWD with the major 
difference being the use of DTIMEL which represeats the user specified 
minimum time-step allowed. The time-step calculations stored in DTIMEU 
proceed exactly as in CNFRWD until the check with DTIMEL is made. If DTIMEU 
(CSGMIN of a node) ^ DTIMEL, the diffusion node temperature calculation is 
identical to CNFRWD. If DTIMEU (CSGMIN of a node) < DTIMEL, the diffusion 
node receives the steady state calculstlon. 

Control constants DTMPCA which contains the allowable diffusion- 
node temperature change and ATMPCA which contains the arithmetic-node 
temperature change are not checked in CNFAST. Thus time-steps are not 
shortened and temperature calculations repeated. The remainder of the 
computational procedure follows those of CNFRWD (Section 6. 3. 1.3). 

6. 3. 2.4 Control Constants 

Control constants DTIMEL, 0UTPUT and TIMEND (>TIME0) must be 
specified as Indicated'in Table 6.2-5 and described in Section 6. 2. 3. 2; 
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otherwise the *'run*' will terminate with an error message. The function of 
optionally specified control constants ASLXCA, BACKUP, DAMPA, DTipiH, NL00P 
and TIME0 is described in Section 6.2. 3.2. As mentioned before in a 
previous paragraph, the user should take considerable amount of caution 
in specifying DTIMEL in order to prevent large pockets of nodes that receive 
the steady state equation Without reiteration. Note also that TIME0 
.may be set negative and that NL00P is set to one if not specified. 

6.3. 2. 5 Error and Other MessaRes 

If control constants DTIMEL, 0UTPUT and TIMEND are not specified, 
the following error message will be printed for each, 

DTIMEL "N0 DTIMEL" 

0UTPUT "N0 0UTPUT INTERVAL" 

TIMEND no message 

A direct check on TIMEND is not made; an indirect message is printed for 
the other explicit routines but is not output for CNEAST. 

If the short pseudo-compute sequence SPCS is not specified, the 
error message will be, 

"CNEAST REQUIRES SHORT PSEUDO-COMPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficient (NDIM < NND) , 
the message will be, 

" LOCATIONS AVAILABLE" 

Note that the nimber printed will be negative indicating the additional 
storage locations required. 

If CSGMIN ^ 0, the message printed will be, 

"C/SK ZER0 or NEGATIVE" 

Checks on the control constants, the pseudo -compute sequence 
and the dynamic storage allocation are made in the following sequence, 
with the rxHi terminating if a single check is not satisfied, 

0UTPUT, DTIMEL, pseudo-compute sequence, and dynamic 

storage locations. 

It should be partictilarly noted that no message is printed if 
ARLXCA is not satisfied with NL00P iterations; ARLXCA and NL00P are 
optionally specified control constants. 
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Table 6.3-2. 


1 . 

2 . 

3. 


4. 

5. 

6 . 

7. 


8 . 


9. 

10 . 

11 . 

12 . 


Basic Computational Steps for GNFAST 

Specification of control constants. Control constants DTIMEL, 0UTPUT 
and TIMEND must be specified. SPCS is required for CNFAST. (Refer to 
Table 6.2-5 for values and Section 6.2, 3. 2 for descriptions.) 

Sufficiency check on dynamic storage. Requirements = NND (NND « 
diffusion nodes) . 

Setting and/or calculation of time-step. At. (Refer to Section 6.2.4 
for detailed procedure.) 

Note that initial time-step equal DTIMEL and subsequent time-step is the 
larger of CSQIIN or DTIMEL. 


Setting of source and diffusion node dynamic storage locations to zero. 
Calling of VARIABLES 1. (Refer to Section 6. 2. 2. 2 for description.) 
Checking of BACKUP. (Refer to Section 6. 2. 3. 2 for description.) 


Calculation of diffusion-node temperatures. (Refer to Section 6. 2. 5.1 
for description.) Calculation differs from the other explicit routines, 
since diffusion nodes with CSGMIN less than DTIMEL receive steady state 
calculation (refer to Section 6. 3. 2. 2.) 

If DTIMEL £ CSGMIN, the node temperature is calculated as, I 

I 


T 

i,n+l 






T. J 


+ 9 . 


i,n 


If DTIMEL > CSGMIN, the node temperature is calculated using the steady 
state expression. 


i,n+l 


q. + Z G.V (T, - T. V 

i»n ij,n j,n i,n^ 


P 

Z 6 
1=1 


ij,n 


Calculation of arithmetic-node temperatures if the number of iterations 
equals NL00P the temperatures are retained without user notification. 
(Refer to Section 6. 2. 5.1 for details.) 

Calling of VARIABLES 2. (Refer to Section 6.2.2. 3 for description.) 


Advancing of time, checking of time to print, and the printing at the 
the output interval. 


Calling of 0UTPUT CALLS. 

Checking for problem end-time stored in user specified control constant 
TIMEND 
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Check and set control constants 


Initialize pointers 

^ ^ ^ 

Initialize and adjust time-step 

T ; ■ : : r ; ■ :;..z 

Set source locations to zero, initialize 
XSP ACE, call VARIALBES 1 


Update Ci, q^, compute Qsum 
and Gsum for each diffusion node 

1 " ' " 

Compute CSGMIN: 

If DTIMEL ^ CSGMIN compute diffusion-node temperature by 
explicit difference method; 

If DTIMEL > CSGMIN compute diffusion-node temperature by 
steady state equation. 

Update and Gj^ once each tim<^step; compute 
arithmetic-node temp, by successive point iteration 

'' "f —————i 

Call VARIABLES 2 
Update time 

I 

If time to print 
call 0UTCAL 



I 



Figure 6.3-4. Functional Flow Chart for CNFAST 
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6.3.3 


Subroutine; CNEXPN 


6. 3.3.1 General Comments 

Subroutine CNEXPN is an explicit routine based upon the exponential 
prediction method;^* the method being imconditionally stable permits any 
size time-steps and requires the short pseudo-compute sequence (SPCS). An 
infinite time-step reduces the transient equation to a steady state one. 
Although the method is unconditionally stable, stability should not be con- 
fused with accuracy. Comparison of several numerical methods, including the 
exponential approximation, is given in Reference 17. 

If accuracy is an important consideration, time-steps should not 
be larger than those taken with the standard explicit method such as used 
in CNFRWD. If high accuracy is not an important consideration, considerable 
savings in computational time can be affected with the use of a large time- 
step. It should be noted that the same savings in computational time may be 
possible with the implicit routines. As another note of interest, CNEXPN 
solutions have a tendency to lag in time the true temperatures. 

6. 3. 3. 2 Finite Difference Approximation and Computational Algorithm 
Diffusion Nodes 

The expression for the numerical method used in subroutine CNEXPN 
for solving the diffusion-node temperatures may he derived from the heat 
balance equation (5.1-6). 


‘^'^i 1 r P 4 4 1 


If G 


ij 


i » 1,2,...,N 
Tj = constant, N < j ^ p 

2 2 

ay + aby (Tj + T^) (Tj + T^) equation (5.1-6) becomes, 


dT 
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i J-f 4. ? 


Gy (I^ - I,) 


1 - 1 , 2... .,11 

Tj = constant, N < j p 


(equation 
5.1-6 of 
Section 5) 


(6.3-6) 
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If we further let G,., q, and T. be invariant with time and 

l3 ^i J 

temperature, equation (6,3-6) may be integrated rather easily to yield, 


-a At 

T *= T e ° + 

i,n+l i,n 


P 

“ .f. ®ij,n ^j,n / 

— ^ \^1 - e “ j (6.3-7) 


2 G.. 

j=l 


where, n = nth time-step; i = 1, 2, N51D (number of diffusion nodes) 

^i*^i’^ij *^ij ~ optionally specified (refer to Tables 6.2-1 - 6.2-4) 

T = constant, (NND + NNA) < j < p (NNA is the number of arithmetic 
nodes and p is the total number of nodes) 

P 

2 G.. 
i=l 



x,n 

At = time-step (refer to Section 6-2.4) 

G,. = a.. +ab.. (T^ + T? )(T. + T. ) 

ij,n ij,n ij,n i,n j,n x,n j,n 


Computationally equation (6.3-7) is applied to the diffusion nodes. 
It should be noted that the foann of equation (6.3-7) represents a "block" 
change in temperatures since the evaluation of n+1 based upon 

Arithmetic Nodes (if any) 

Arithmetic-node temperatures are calculated in the same manner as 
in CNFRWD (Section 6. 3. 1.2) or refer to Section 5.2.3 for finite difference 
algorithm. 


6.3. 3. 3 Comments on the Computational Procedure 

The important steps of the computational procedure used in sub- 
routine CNEXPN are indicated in Table 6.3-3 and a functional flow chart is 
shown in Figure 6.3-5. A detailed computational procedure requires the 
examination of the CNEXPN computer listing which is presented in Appendix A 
but some general computational details are given in Section 6. 2. 5.1. The 
computational process of subroutine CNEXPN is essentially identical to 
CNFRWD with the difference being the finite difference expression used for 
the calculation of the diffusion nodes and the time-step which is calculated 
as CSGMIN*CSGFAC in lieu of CSQIIN/CSGFAC. The "look ahead" feature for 
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tiaie-step calculation as well as a check with DTIMEH, DTIMEL and DIMPCA is 
identical to CNFB.W). T^peratures of arithmetic nodes are calculated after 
the diffusion nodes and utilize M400P, AKLXCA, and DAMOPA in exactly the 
same way as CNFRWD. The verbal flow description of GKFRWD (Section 6. 3. 1.3) 
applies here except for the differences indicated above, 

6. 3. 3.4 Control Constants 

Control constants 0UTPUT and TIMEND (> TIME0) must be specified 
as Indicated in Table 6.2-5 and described in Section 6.2. 3.2; otherwise 
the "run" will terminate with an error message. The function of optionally 
specified control constants APiXCA, ATMPCA, BACKUP, CSGFAC, DAMPA, DTIMEH, 
DTIMEL, DTMPCA, NL00P, and TIME0 is described in Section 6. 2. 3. 2. The user 
should take particular care in the selection of CSGFAC since too large of 
a time-step would lead to grossly inaccurate temperatures even though the 
solution is stable. Note also that TIME0 may be set negative and that 
NL00P is set to one if not specified. 

6,3,3. 3 Error and Other Messages 

If control constants 0DTPUT and TIMEND are not specified, the 
following error message will be printed for each, 

0UTPUT "N0 0UTPUT INTERVAL" 

TIMEND "TIME STEP T00 SMALL" 

The reason for the TIMEND error message is that a direct check on TIMEND 
is not made; the resultant error message just happens to be a quirk in the 
coding. 

If the short pseudo-compute sequence SPCS is not specified, the 
error message will be, 

"CNE3CPN REQUIRES SHORT PSEUDO-COMPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficient 
NDIM < (NND + NNA)), the message will be, 

" LOCATIONS AVAILABLE" 

Note that the number printed will be negative Indicating the additional 
storage locations required. 
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If the time-step used is less than the time-step allowed 
(DTIMEL) which may be optionally specified by the user, the message will be, 

"TDIE STEP T00 SMALL" 

If CSGMIN < 0, the message printed will be, 

"CSGMIN ZER0 or NEGATIVE" 

Checks on the control constants, the pseudo-compute sequence and 
the dynamic storage allocation are made in the following sequence with th® 
run terminating if a single check is not satisfied, 

0UTPUT, pseudo-compute sequence, dynamic storage locations. 

It should be particularly noted that no message is printed if 
ASLXCA is not satisfied with NL00P iterations; AELXCA and NL00P are 
optionally specified control constants. 



Table 6.3-3. 



2 . 

3. 

4. 

5. 

6 . 
7, 




8 . 

9. 

10 . 

11 . 

12 . 

13. 

14. 


Basic Coii^utational Steps for GNEXPN 

Specification of control constants. Control constants 0UTPUT and TIMEND 
must be specified. SPCS is required for GNEXPN. (Refer to Table 6.2-5 
for values and Section 6. 2. 3. 2 for description.) 

Sufficiency check pn dynamic storage. Requirements = NND + NNA (NND = 
diffusion nodes and NNA = arithmetic nodes). 

Setting and/or calculation of time-step, At. (Refer to Section 6.2,4 
for detailed procedure,) Calculated time-step = 0.95 * CSGMIN * CSGFAC. 

Setting of source and diffusion node dynamic storage locations to zero. 

Calling of VARIABLES 1. (Refer to Section 6. 2. 2. 2 for description.) 

Checking of BACEUP. (Refer to Section 6 . 2 . 3. 2 for description.) 

Calculation of diffusion-node temperatures. (Refer to Section 6. 2, 5.1 
for description.) 

Diffusion-node temperatures are calculated by using (refer to 
Section 6. 3. 3, 2), 


-“ At 'l,n ''' ‘"Id.” 

'l.n+1 * \,n ^ 
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Erasure of all temperature calculations for latest time-step if allowable 
temperature change criterion DTMPCA is not satisfied and recalculation 
of temperatures with reduced time-step. 


Calculation of arithmetic-node temperatures. If the number of iterations 
equal NL00P, the temperatures are retained without user notification 

Erasure of all temperature calculations for latest time-step if allowable 
temperature change criterion ATMPCA is not satisfied and recalculation 
of temperatures with reduced time-step. 

Calling of VARIABLES 2 and checking of BACKUP. (Refer to Section 6. 2. 2. 3 
and 6. 2. 3. 2 for description.) 

Advancing of time, checking of time to print, and the printing at the 
output interval. 

Calling of 0UTPUT CALLS. . 

Checking for problem end time stored in user specified control constant 
TIMEND. 
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6.3.4 


Subroutine; CNDUFR 


' 6. 3. 4.1 General Comments 


Subroutine CNDUFR is an explicit numerical solution routine that 
uses an tmconditionally stable DuFort-Frankel method.^* 12,17 
DuFort-Frankel method replaces the present tenperature of the node being 
operated on by the average of future and past temperatures in the forward 
differencing equation. In subroutine CNDUFR the present temperature of the 
node being operated on is replaced by a time-weighted average of future and 
past temperatures , CNDUFR requires the short pseudo-compute sequence (SPCS) , 

The intent of an unconditionally stable routine such as CNDUFR 
is the reduction of computational time by using time-steps greater than 
those allowed with the conditionally stable explicit methods as constrained 
by the stability criterion. However, less accuracy can be expected with a 
lengthened time-step. The time-step controlled with control constant CSGFAC 
. represents a user decision that is difficult and must be aided by a trial 
and error procedure. 


Examination of several CNDUFR solutions reveals a tendency to 
) lead in time the true temperatures. 


6. 3. 4. 2 Finite Difference Approximation and Computational Algorithm 

Diffusion Nodes 

The DuFort-Frankel explicit finite difference expression®* 
for calculating the diffusion-node temperatures may be readily determined 
as follows ; 


Using the standard explicit finite difference expression. 


^^i,n+l ~ ^i.n ^ 
At 
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q. + 2 G.. (T. - T. ) 

^i,n X 3 ,n 3 ,n x,n' 


(6.3-8) 


i « 1,2,. ..,N 
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where, ** i =* 1,2, ...,M (eqiial titee-steps) 

and defining. 


“ C^/At (refer to Section 6.2.4 for discussion on At) (6.3-10) 
equation (6. 3-^8) can be expressed as. 


C. T. , + 2 q. + Z G. . (2 T. - T. ,) 

i,n i,n-l ^i,n i 3 »n 3 ,n x,n-l' 


i,n+l 


(6.3-11) 


^ij ,n 


i * 1,2,..., N 


In CNDUFR the present temperature, T, , of equation (6.3-8) is replaced 

x,n 

by a weighted average of future temperature, n+1* past temperature, 
Tf The weighting is based bn unequal time-steps. 

(At , T, + At T, -) 

T * n-1 x,n+l n x,n-l 

^i,n ” At 1 + At 

n-1 n 

where. At , = t - t , (past time-step) 
n— X n n-1 


(6.3-12) 
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(6.3-13) 

(6.3-14) 


Equation (6.3-8) becomes, - 
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(6.3-15) 


where, i =* 1,2,..., NND (number of diffusion nodes) 

T. = constant, (NND + NNA) < J < p (^A is the number of arithmetic 

j 

nodes and p is the total number of nodes) 

G.. = a.. +ab.. (T? + T? )(T. + T . ) 

ij,n xj,n xj ,n ' 3 ,n x,n' ' j ,n i,n'' 

^i’^^i’^ij’^ij * optionally specified (refer to Tables 6.2-1 - 6.2-4) 
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In CNDUFR, equation (6 . 3-15) is applied to the diffusion nodes 
with the computational procedure being a "block" change in temperature from 
one time-step to another. 

Arithmetic Nodes (if any) 

Arithmetic-node temperatures are calculated in the same manner 
as in CNFRTO (Section 6. 3. 1.2) or refer to Section 5.2.3 for the finite 
difference algorithm. 

6. 3. 4. 3 Comments on the Computational Procedure 

The important steps of the computation procedure used in sub- 
routine CNDUFR are indicated in Table 6.3-4 and a functional flow chart 
is shown in Figure 6.3-6. A computer listing of CNDUFR is fotand in 
Appendix A but some general computational details are given in Section 6. 2. 5.1. 
Tlie computational procedure for CNDUFR follows the CNEXPN computational 
pattern, but with the temperatures of the diffusion nodes calculated by 
the DuFort-Frankel method of the exponential prediction method. Another 
significant difference is that CNDUFR must provide for two sets of past 
temperatures which are required for DuFort-Frankel algorithm; two time- 
steps for consecutive time-step calculations are also required. Otherwise, 
checks and control constant use are identical to CNEXPN. Thus, the verbal 
flow description of Section 6. 3.1. 3 applies directly except for the 
differences indicated above. 

6. 3.4.4 Control Constants 

Control constants 0UTPUT and TIMEND (> TIME0) must be 
specified as indicated in Table 6.2-5 and described in Section 6. 2. 3. 2; 
cbhsrwise the run will terminate with an error message. The function 
of optionally specified control constants AELXCA, ATMPCA, BACKUP, CSGFAC, 

DAMPA, DTIMEH, DTIMEL, DTMPCA, NL00P, and TIME0 is described in 
Section 6.2. 3.2. The user should take particular care in the selection 
of CSGFAC since too large. of a time-step would lead to grossly inaccurate 
temperatures even though the solution is stable. Note also that TIME0 
may be set negative and that NL00P is set to one if not specified. 
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6. 3. 4. 5 Error and Other Messages 


If control constants 0UTPUT and TIMENI) are not specified, the 
following error inessage will be printed for each, 

0UTPUT *'N0 0UTPUT INTERVAL" 

TIMEND "TIME STEP T00 SMALL" 

The reason for the TIMEND error message is that a direct check on TIMEND 
is not made; the resultant error message just happens to be a quirk in 
the coding. 

If the short pseudo-compute sequence SPCS is not specified, the 
error message will bej < 

"CNDUFR REQUIRES SHORT PSEUDO-COMPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficient 
(NDIM < (2*NND + NNA)), the message will be, 

LOCATIONS AVAILABLE" 

Note that the number printed will be negative indicating the additional 
storage locations required. 

If the time-step used is less than the time-step allowed (DTIMEL) , 
which may be optionally specified by the user, the message will be, 

"TIME STEP T00 SMALL" 

If CSGMIN £ 0, the message printed will be, 

"CSGMIN ZER0 or NEGATIVE" 

Checks on the control constants, the pseudo-compute sequence and 
the dynamic storage allocation are made in the following sequence with the 
run terminating if a single check is not satisfied, 

0DTPUT, pseudo-compute sequence, dynamic storage locations 

It should be particularly noted that no message is printed if 
ARLXCA is not satisfied with NL00P iterations; AELXCA and NL00P are 
optionally specified control constants. 
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Table 6.3-*4 


1 . 

2 . 

3 . 

4. 

5 . 

6 . 
7 . 


8 . 

9 , 

10 . 

11 . 

12 , 

13. 

14. 


Basic Computational Steps for CNDUFR 

Setting of control constants to nominal values. Control constants 
0DTPUT and TIMEND must be specified. SPCS is required for CNDUFR. 
(Refer to Table 6.2-5 for values and Section 6. 2. 3. 2 for description.) 

Sufficiency check on d 3 mamic storage. Requirements = 2*NND + HNA 
(NND = diffusion nodes and MA = arithmetic nodes). 

Setting and/or calculation of time-step. At. (Refer to Section 6.2.4 
for detailed procedure.) Calculated time-step = 0.95*CSGMIN*CSGFAC. 

Setting of source and diffusion node dynamic storage locations to zero. 

Calling of VARIABLES 1. (Refer to Section 6. 2. 2. 2 for description.) 

Checking of BACKUP. (Refer to Section 6.2. 3.2 for description.) 

Calculation of diffusion-node temperatures. (Refer to Section 6. 2, 5.1 
for description.) 

Diffusion-node temperatures are calculated by using (refer to 
Section 6. 3. 4. 2), 


_ P P 

T T. i (C. - I G.. ) + 2 G.. T. + q. 

n i,n-l i,n X 3 ,n xj »n j ,n ^i,n 

'‘i,n+l _ _ p 

C - T , (C. - 2 G.. ) 

l,n n-1 ' x,n xj,n' 


where. 
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h-1 


n-1 At , + At 
n-1 n 


At 

- = P 

n At 1 + At 
n-1 n 

Erasure of all temperature calculations for latest time-step if 
allowable temperature change criterion DTMPCA is not satisfied and 
temperature recalculation with reduced time-step. 

Calculation of arithmetic-node temperatures; if the number of itera- 
tions equal NL00P , the temperatures are retained without user 
notification (refer to Section 6.2. 5.1 for details). 

Erasure of arithmetic-node temperatures for latest time-step if allowable 
temperature change criterion A1MPCA is not satisfied and temperature 
recalculation with reduced time-step. 

Calling of VARIABLES 2 and checking of BACKUP. (Refer to Section 6. 2. 2. 3 
and 6. 2. 3. 2 for description.) 

Advancing of time, checking of time to print, and the printing at the 
output interval. 

Calling of 0UTPUT CALLS. 

Checking for problem end time stored in user specified control constant 
TIMEND. 
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Figure 6.3-6. 


Functional Flow Chart for CNDUFR 
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6-3.5 


Subroutine; CNQUIK 
6. 3. 5.1 General Comments 

Subroutine CNQUIK is a numerical solution routine that uses an 
algorithm composed of half DuFort-Frakel method®* and half 

eiq>onentlal prediction method.^* CNQUIK requires the short pseudo- 
compute sequence (SPCS)j characteristics of subroutines CNDUiK and CNEXPN, 
as described in Section 6.3.3 and 6.3.4, also apply to CNQUIK. 

Why CNQUIK? Examination of CNDUFR and CNEXPN solutions reveals 
that CNDUFR has a tendency to yield temperatures which lead the true tempera- 
tures, whereas CNEXPN has a tendency to lag the true temperatures. Thus, 
it was theorized that a combination of CNDUFR and CNEXPN should yield a 
more accurate solution than either one. Preliminary results Indicate that 
CNQUIK is more accurate than either CNDUFR or CNEXPN with approximately 
the same solution time. It can also be theorized that a more accurate 
combination of the DuFort-Frankel and exponential prediction is probably 
possible than the half and half used in CNQUIK. However, a detailed 
study will be required before a realistic evaluation of CNQUIK can be made. 


6.3. 5. 2 Finite Difference Approximation and Computational Algorithm 

Diffusion Nodes 


Subroutine CNQUIK uses a numerical solution algorithm composed 
of half DuFort-Frankel and half exponential prediction. That is the 
temperature of the diffusion nodes is calculated by using. 
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n « nth time-step 

i * 1,2,...,NHD (number of diffusion nodes) 

T. “ constant, (NND + NNA) < j < p (NNA is the ntnnber of 
3 

arithmetic nodes and p is the total number of nodes) 


a 


P 

Z G 


ij,n 


"i,n 


Gj. = a.. + ab.. (T? + T? )(T. + T, ) 

ij,n l3,n ij,n ■ i,n j»n j»a 


At 


At 


h-1 


n At , + At ’ n-1 At , + At 

n-1 n n-1 n 


C^, q^, “ optionally specified (refer to Tables 6.2-1 - 6.2-4) 


C. = C^/At (refer to Section 6.2.4 for discussion of At) 

Arithmetic Nodes (if any) 

Temperatures of arithmetic nodes are calculated in the same manner 
as in CNFRWD (Section 6.3. 1.2) or refer to Section 5.2.3 for the finite 
difference algorithm. 

6.3.5. 3 Comments on the Computational Procedure 

The important steps of the computational procedure used in sub- 
routine CNQUIK are Indicated in Table 6.3-5 and a functional flow chart is 
shown iu Figure 6.3-7. A computer listing of CNQUIK is foxmd in Appendix A. 
General computational details are given in Section 6.2. The computational 
procedure for CNQUIK follows CNEXPN or CNDUFE with the diffusion-node 
temperatures Calculated with the half DuFort-Frankel and half exponential 
prediction algorithm being the only difference. Arithmetic-node tempera- 
tures are calculated in the same manner as the other SINDA explicit 
routines. Note that the time-step is calculated as CSGMIN*CSGFAC and checks 
are the same as CNEXPN or CNDUFR. Thus, the verbal flow description of 
Section 6. 3.1. 3 applies directly except for the differences indicated 
above. 
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^ 6. 3. 5. 4 Coatrol Constants 

Control constants 0UTPUT and TIMEND (> TIME0) must be specified 
as indicated in Table 6.2-5 and described in Section 6. 2. 3. 2; otherwise 
the "run” will terminate with an error message. The function of optionally 
specified control constants ARLXCA, ATMPCA, BACiaJP, CSGFAC, DAMPA, DTiMEH, 
DTIMEL, DTMPCA, NL00P, and TIME0 is described in Section 6. 2. 3. 2. Again, 
caution must be exercised in the selection of CSGFAC since too large of a 
time-step would lead to grossly inaccurate temperatures even though the 
solution is stable. Note also that TIME0 may be set negative and that 
NL00P is set to one if not specified. 

6. 3. 5. 5 Error and Other Messages 

If control constants 0UTPUT and TIMEND are not specified, the 
following error message will be printed for each, 

0UTPUT "N0 0DTPUT INTERVAL" 

TIMEND "TIME STEP T00 SMALL" 

The reason for the TIMEND error message is that a direct check on TIMEND 
is not madei the resultant error message just happens to be a quirk in 
the coding. 

If the short pseudo-compute SPCS is not specified, the error 
message will be, 

"CNQUIK REQUIRES SHORT PSEUDO-COMPUTE SEQUENCE" 

If the d 3 mamic storage allocation is not sufficient 
(NDIH < (2*NND + NNA), the message will be, 

LOCATIONS AVAILABLE" 

Note that the number printed will be negative indicating the additional 
storage locations required. 

If the time-step used is less than the time-step allowed (DTIMEL) , 
which may be optionally specified by the user, the message will be, 

"TIME STEP T00 SMALL" 

If CSGMIN ^ 0, the message printed will be, 

"CSGMIN ZER0 or NEGATIVE" 
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Checks on the control constants, the pseudo-compute sequence 
and the dynamic storage allocation are made in the following sequence with 
the inxn terminating if a single check is not satisfied, 

0UTPUT, pseudo-compute sequence, dynamic storage locations. 

It should be particularly noted that no message is printed if 
ASLXCA is not satisfied with NL00P iterations; AELXCA and NL00P are 
optionally specified control constants. 
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Table 6.3-5. Basic Computational Steps for CNQUIK 

1. Specification of control constants. Control constants 0UTPUT and 
TIMENB must be specified. SPCS is required for CNEXPN. (Refer to 
Table 6.2-5 for values and Section 6.2. 3.2 for description.) 

2. Sufficiency check on d 3 mamic storage. Requirements = 2(NND) + NNA 

* diffusion nodes and NKA. = aritlimetic nodes). 

3. Setting and/or calculation of time-step, At. (Refer to Section 6.2.4 
for detailed procedure.) Calculated time-step = 0.95*CSGMIN*CSGFAC. 

4. Setting of source and diffusion node dynamic storage locations to zero. 

5. Calling of VARIABLES 1. (Refer to Section 6. 2. 2. 2 for description.) 

6. Checking of BACKUP. (Refer to Section 6. 2. 3. 2 for description.) 

7. Calculation of diffusion-node temperatures. (Refer to Section 6. 2. 5.1 
for description.) 

Diffusion-node temperatures are calculated by using (refer to 
Section 6. 3.5.2), 

^i,n+l “ ^^CKDUFR ^ ^CNEXPN^^^‘° 

(Refer to equation 6.3-17, Section 6. 3. 5. 2.) 

8. Erasure of all temperature calculations for latest time-step if 
allowable temperature change criterion DTMPCA is not satisfied and 
temperature recalculation with reduced time-step. 

9. Calculation of arithmetic-node temperatures; if the number of 
iterations equal NL00P, the temperatures are retained without user 
notification. (Refer to Section 6. 2. 5.1 for details.) 

10. Erasure of all temperature calculations for latest time-step if allow- 
able temperature change criterion AIMPCA is not satisfied and tempera- 
ture recalculation with reduced time-step. 

11. Calling of VARIABLES 2 and checking of BACKUP. (Refer to Section 
6. 2. 2. 3 and 6. 2. 3. 2 for description.) 

12. Advancing of time, checking of time to print, and the printing at the 
output interval. 

13. Calling of 0DTPUT CALLS. 

14. Checking for problem end time stored in user specified control constant 
TIMEND. 
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Check an«J s-et control constants 


Initialize pointers 


Initialize and adjust time-step 


Set source locations to zero; initialize 
XSPACE, call VARIABLES 1 


Update Gi, q^, G^; compute and 

Ggujjj for each diffusion node 


Compute CSQIIN; compute diffusion-node temperatures 
by half DuFort-Frankel and half exponential 
prediction method; compute DTMPCC 



1st pass 


3TMPCC olc 


Check and adjust 
time-step 


Undo diffusion-node I ^ iShorten 

temperature calculations I ' [time-step 



Set time-step 


Update qi & once each time-step. 
Compute arithmetic-node temperatures 
by successive point iteration 


Compute ATMPCC 


DTMPCC oiT 


BACKUP 0^ 


Call VARIABLES 2 


flMEND 


If time to print 
call 0UTCAL 


Update time 


Return 


Figure 6.3-7. Functional Flow Chart for CNQUIK 
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6.4 


Transient Implicit Solution Soutines 

SINDA implicit solution routines number three; these routines 
are Identified as follows: 

CNBACK Implicit backward difference method. 

Requires long pseudo-compute sequence (LPCS) . 

CHFWBK Implicit forward-backward differencing, using 

Crank-Nicolson method. 

Requires long pseudo-compute sequence (LPCS). 

CNVARB Combination of CNBACK and CNFWBK. 

Requires long pseudo-compute sequence (LPCS), 

Implicit methods generally tend to be more accurate than explicit 
methods and are unconditionally stable as are some explicit methods. With 
Implicit methods the time-step is specified in contrast to the calculated 
time-steps of explicit methods with their stability criterion. An important 
consideration in the use of implicit methods is that the time-step DTIMEI 
should be specified in conjunction with control constant NL00P which 
represents the maximum number of computational iterations during each time- 
step. Since each iterative calculation is essentially equivalent to a 
time-step calculation for an explicit method, the combination of DTIMEI and 
ML00P for a given time period should be set less than the total number of 
time-steps used by the explicit method for the same time period, if com- 
putational time is to be reduced; this of course assumes that during each 
time-step the maximum number of iterations is required. If the NL00P 
iterations are required during a time-step, the temperature accuracy is 
affected but the magnitude would depend upon the value used for the maximum 
allowable relaxation temperature change criteria, ARLXCA and DRLXCA. It 
should be noted if NL00F iterations are required during a time-step, the 
message "RELAXATION CRITERIA NOT MET" is printed. 

A detailed description of each implicit routine, as presented 
on the pages to follow, relies on the general description of Section 6.2. A 
brief description of these routines is summarized first . 

CNBACK Uses the standard backward differencing algoritlun and 
requires the long pseudo-compute sequence (LPCS) . The time-step must be 
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specified via control constant DTIMEI and used in conjxmction with the con- 
trol constant NL00P. CNBACK uses the acceleration of convergence feature. 

CNFHBK uses the Cranfc-Nicolson algorithm which is composed of 
half forward differencing and half backward differencing. CNFWBK solutions 
tend to be more accurate than CNMCK solutions with approximately 25% less 
iterations; however CNFWBK solutions have "blown” on occasions. 

CNVARB uses a combination of forward differencing and backward 
differencing. Unlike CNFWBK which is half and half, CNVAKB uses a variable 
beta factor which ranges from 0 to 1 . Thus CNVARB uses a method that is 
somewhere between forward differencing and backward differencing. 
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6.4.x 


Subroutine; CNBACK 


6. 4. 1.1 General Gotaments 

Subroutine CNBACK is an implicit routine that uses the standard 
backward difference expression and requires the long pseudo-compute (LPCS) . 
Time-step must be specified via DTIMEI otherwise the "run" will terminate 
with an error message printout. The time-step value is arbitrary but the 
user should consider DTHiEI in conjunction with the control constant NL00P 
which represents the maximum number of ccmiputational iterations during each 
time-step (refer to Section 6. 2. 3. 2 for description). 

Implicit methods tend to be more accurate than explicit methods 
and are unconditionally stable, but Implicit solutions often oscillate at 
start up or boundary step changes when heat transfer by radiation is 
present. CNBACK internally controls sudden radiation beat transfer 
changes by an averaging technique which is termed, "radiation damping” (refer 
to Section 6.2.6 for details). This automatic damping has been very effec- 
tive in many solutions that have been examined and lessens the need for the 
use of DAMPD and DAMPA. 

6. 4. 1.2 Finite Difference Approximation and Computational Algorithm 

The numerical solution algorithm used in subroutine CNBACK is 
the standard backward-difference expression^^ * which may be expressed 

as: 


(T 


i.n-fl 


l.n 


At 


% t, 

i,n 
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j*l 
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ij ^^j.n+1 


■ ^i,n+l^ 
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j»n+l 


i,n+l 


) 


(equation 5.2-5 of Section 5.2.2) 
i » 1,2, . . . ,N 

r. ,, * constant, N < j < p 
3,n-ri — 


The computational procedure for the backward difference formula- 
tion vmst necessarily be re-iterative because of the need to solve a set of 
simultaneous nonlinear equations. 
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Diffusion Nodes 


Diffusion node temperatures are solved by "successive point" 
iteration but differs from the arithmetic-node temperature calculation 
because of the capacitance term and the use of "radiation damping" (refer 
to Section 6. 2. 5. 2). 


T 

i,k+l 


DD* T 


i,k 


+ DN* 


_ i P 

C. T. + q. + 2 6.. T. , + 2 G.. T. , - (qj 

i,n i,n ^i,n ij,n j,k+l . ij ,n j ,k ^i'ave 

■ ■ ■ - ■ ■■ ■ - ■ . ■ . :■ . ■ 


p 

+ 2 a^ 


(6,4- 


i»n ij.n 


where, i = 1,2,... , NND 

n = nth time-step 
k “ kth iteration 

DN * DAMPD (diffusion-node damping factor) 

DD = 1.0 - DN 

3 

®ij ,n “ *ij ,n ^j ,A (il “ k if j > i and I = k+1 if j < i) 

^i**^i*^ij "'^ij * optionally specified (refer to Tables 6,2-1 - 6.2-4) 

C. =» C. /At (At « time step, refer to Section 6.2.4) 

P 4 4, 

(qi)ave “ ^ij,n *-^^i,k^ ^ (T^^j^) 2 ]/ 2 . 0 , average heat loss from 

the 1th node (refer to Section 6.2.6 on radiation damping for 
details) 


Details on the computational procedure for implicit routines are 
presented in Sections 5.2.2 and 6.2. 5. 2. 

Arithmetic Nodes 

Arithmetic-node temperatures are calculated Identically the same 
in all the SINDA numerical solution routines. Thus, refer to either 
Section 6.3. 1.2 or Section 6.2. 5. 2 for the finite difference algorithm. 

6.4. 1.3 Comments on the Computational Procedure 

The important steps of the computational procedure used in sub- 
routine CNBACK are indicated in Table 6.4-1. For a detailed step-by-step 
computational description, the user must examine the computer listing for 
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CNBACK In Appendix B, but some general computational details are given in 
Section 6. 2. 5. 2. A functional flow chart of CNBACK is shown in Figure 6.4-1. 

Both diffusion-node tanperatures and arithmetic-node temperatures 
arc calculated by "successive point" iteration. Each third iteration, 
diffusion-node temperatures which are decreasing over two time-steps are 
extrapolated in an attempt to accelerate convergence (refer to Section 6.2.7). 
Temperature convergence is examined during each time-step by checking DRLXCG 
and AEJJCCC against the user control constants DELXCA (for diffusion nodes) 
and ARLXCA (for arithmetic nodes), respectively. If temperatures have not 
converged with NL00P iterations, the message "SEIAXATION CRITERIA NOT MET" 
is printed. Control constant NL00P is used to specify the maximum nximber of 
it nations allowed during each time-step. 

VARIABLES 1 and VARIABLES 2 are performed only once for each time- 
step. Since this subroutine is implicit, the user must specify the time- 
step to be used through the control constant DTIMEI in addition to control 
constant TIMEND and 0UTPUT. The look ahead feature for the time-step 
calculation used in CNFRWD is also employed in CNBACK as are checks for 
maximum allowable time-step DTIMEH, maximum allowable temperature change 
between time-steps^ DTMPCA (diffusion nodes) and ATMPCA (arithmetic nodes). 
The minimum time-step DTIMEL is not checked however. Damping of solutions 
can be achieved through the use of the control constants DAMPD and DAMPA 
but "radiation damping" (refer to Section 6.2.6) used by CNBACK lessens the 
need for the damping factors DAMPD and DAMPA. 

6. 4.1. 4 Control Constants 

Control constants ARLXCA, DRLXGA, DTIMEI, NL00P, 0UTPUT, and 
TIMEND must be specified as indicated in Table 6.2-5 and as described in 
Section 6. 2. 3. 2; otherwise "run" will terminate with an appropriate error 
message. The function of optionally specified control constants ATMPCA, 
BACKUP, DAMPA, DAMPD, DTIMEH, DTMPCA, and TIME0 is described in 
Section 6. 2. 3. 2. 

Specification of time-step DTIMEI should be done in conjunction 
with control constant NL00P which represents the maximum number of com- 
putational iterations during each time-step. Since each iterative calcula- 
tion is essentially equivalent to a time-step calculation for an explicit 
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method, the combination of DTBffil and NL00P for a given time period should 
be less than the total number of time~steps by the explicit method for the 
same period. Note also that TIME0 may be set negative. Specification of 
ARLXCA and DELXCA depends upon the problem but a typical value is 0.1. 

6.4.1. 5 Error and Other Messages 

If control constants MLXCA, DKIiXCA, DTIMEI, NL00P, 0UTPUT and 
TIMEND are not specified, the following error message will be printed for 
each. 


ARLXCA 

"N0 ARLXCA” 

DRLXCA 

”N0 DRLXCA" 

DTIMEI 

”N0 DTIMEI” 

NL00P 

”N0 NL00P” 

0UTPUT 

”N0 0UTPUT INTERVAL” 

TIMEND 

"TRANSIENT TIME N0T SPECIFIED' 


If the long pseudo-compute sequence LPCS is not specified, the 
error message will be, 

"CNMCK REQUIRES LONG PSEUDO-COMPUTE SEQUENCE” 

If the dynamic storage allocation is not sufficient 
(NDIM < (3*NND + NNA + NNB)), the message will be, 

” LOCATIONS AVAILAH.E" 

Note that the ntimber printed will be negative indicating the additional 
storage locations required. 

If CSGMIN ^ 0, the following message will be printed, 

"CSrailN ZER0 or NEGATIVE” 

If either ARLXCA or DRI^CA is not satisfied with NL00P iterations, 
the following message will be printed, 

"RELAXATI0N CRITERIA N0T MET” 

Checks on the control constants, the pseudo-compute sequence 
and the dynamic storage allocation are made in the following sequence with 
the "run” terminating if a single check is not satisfied, 

NL00P, TIMEND, 0UTPUT, ARLXCA, DTIMEI, DRLXCA, LPCS and dynamic 
storage allocation. 
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Table 6.4-1. Basic Computational Steps for CNBACK 

1. Specification of control constants. Control constants ARLXCA (if 
NNA > 0), DRLXCA (if NOT) > 0), DTIMEI, NI00P, 0UTPUT and TIHEND 
(TIMESD _> TIME0) must be specified. LPCS is required. (Refer to 
Table 6.2-5 for values and Section 6. 2. 3. 2 for description.) 

2. Sufficiency check on dynamic storage. Requirements = 3*NND + NNA + NNB 
(NOT) * diffusion nodes, NNA = arithmetic nodes and NNB = bowdary nodes). 

3. Setting and/or calculation of time-step, At. (Refer to Section 6.2.4 
for detailed procedure.) Time-step = DTIMEI. 

4. Setting of iterative D0 loop, 1 to NL00P. 

5. Setting of source locations to zero. 

6. Calling of Variables 1. (Refer to Section 6. 2. 2, 2 for description.) 

7. Checking of BACKUP. (Refer to Section 6. 2. 3. 2 for description.) 

8. Diffusion-node temperature calculations, first Iteration only. 

Evaluation of q^, and 

Damping of radiation heat transfer. (Refer to Section 6. 2. 5. 2.) 

Calculation of diffusion-node tenperature. 

The computational algorithm depends upon the presence of radiation 

heat transfer, but the method of solution is the standard implicit 

algorithm (refer to Section 6. 2. 5. 2). 

9. Conversion of T^ to degrees Rankine. 

10. Diffusion-node tenperature calculations, successive iterations after first. 

Repeating of step 8, except that q^^ C^ and are not updated. 

Calculation of DRIiXCC. 

11. Acceleration of convergence every third iteration if linear extrapola- 
tion is met (refer to Section 6.2.7). 

12. Conversion of T^ to degrees Fahrenheit* 

13. Calculation of arithmetic-node temperatures, second and succeeding 
Iterations; arithmetic-node temperatures are not calculated on the 
first Iteration (refer to Section 6. 2. 5. 2 for details). 

14. Conversion of temperatures to degrees Rankine. 

15* Checking of ARLXCA and DRLXCA for convergence and 0PEITR for output. 

If both ARLXCA and DRLXCA are satisfied, iterations during a time-step 
ceases, otherwise NL00P iterations are performed. 

16. Checking of ATMPCA and DTMPCA* If either one is not satisf ied time-step 
is shortened, previous. temperatures erased, and temperatures recalculated 
for shortened time-steps (refer to Section 6. 2. 5. 2). 

17. Conversion of temperatures back to degrees Fahrenheit. 

18. Calling of VARIABLES 2 and checking of BACKUP (refer to Section 6. 2. 2. 3 
and 6. 2. 3.2). 

19. Advancing of time, checking of time to print, and the printing of the 
output interval. 

20. Calling of 0UTPUT CALLS. 

21. Checking for problem end time stored in control constant TIMEND. 
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6 . 4.2 Subroutine; CNFWBK 

6.4.2. 1 General Conments 

Subroutine CNFWBK is an Implicit numerical solution routine that 
uses the Crank-Nocolson algorithm.^* ®» The long pseudo-compute sequence 
(LPCS) is required and the nodal temperatures (both diffusion and arithmetic) 
are solved by "successive point" iterations. The iteration looping, con- 
vergence criteria and other control constant checks are identical to CNBACK, 
Time-step must be specified via control constant DTlMEI. Diffusion and 
arithmetic temperature calculations may be damped through the use of DAMPD 
and DAMPA, respectively. Thermal radiation heat transfer is uniquely 
"handled" via a so-called "radiation damping" (refer to Section 6.2.6), 
and acceleration of convergence (refer to Section 6.2.7) is also available 
la CNFWBK. 


CNFWBK solutions which are based on a half forward differencing 
and a half backward differencing method tend to be more accurate than CNBACK 
solutions with approximately the same solution time. 


6. 4. 2. 2 Finite Difference Approximation and Computational Algorithm 


The numerical solution algorithm used in subroutine CNFWBK is the 
Crank-Nicolson method, which is half forward differencing and half backward 
differencing, and may be expressed as: 


CT. - T. ) , 

- l,n-<-l x,n _ 1 . _ . 

^i At 2 '■‘■forward ■‘‘backward'^ 


(6.4-2) 


^forward - ’l,n HurPi.n "l,n> * 


■'tacfcuard " ’l,n ®lj,n '^i.a+1^ “’’ij.n ^^3,n+l" 

n » nth time-step 

1 » 1,2,...,N 

p =» total number of nodes 

■'j.n' ^j,n+l * N < 3 < p 
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The computational procedure for the forward-backward differaice 
formulation must be re-iterative because of the need to solve a set of 
simultaneous nonlinear equations. The pattern of computation is very 
similar to that used in CNBACK. 


Diffusion-Nodes 

Diffusion node temperatures are solved by ’'successive point" 
iteration but the algorithm differs from the algorithm used in CNBACK 
because of the additional terms arising from the forward difference 
portion of the expression. 


where. 


T 

l,k+l 

» DD* T. . + DN* [Q - (q.) ]/G 

l,k stom 1 ave sum 

(6.4-5) 
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Qi = 2 


- P 

2 G. T. + Z a.. (T. - T. ) 

x,n x,n xj,n j,n i,n' 

(6.4-7) 

G = 

sum 


P 

+ Z a . j 
j-i 

(6.4-8) 


n * nth time-step; k *= kth iteration 
^i’*^i*^ij * optionally specified (refer to Tables 6.2-1 - 6.2-4) 

DN DAMPD (diffusion-node damping factor) 

DD « 1.0 - DN 

C. * C, /At (At “ time-step) 

^‘^i^ave “ ^ ^^ij,n *‘^'*^i,k^ ^ heat loss from 

ith node (refer to Section 6.2.6 on radiation damping 
for details) 

(Note that the known quantities at time-step, n, are indicated by Q^, 
equation 6.4-7.) 
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Arithmetic Nodes 


Arithmetic-node temperatures are calculated identically the same 
in all the SINDA numerical solution routines. Thus, refer to Section 

6.3.1.2 or Section 6. 2.5. 2 for the finite difference algorithm. 

6. 4.2. 3 Comments on the Computational Procedure 

The important steps of the computational procedure used in sub- 
routine CNFWBK are indicated in Table 6.4-2, For a detailed step-by-step 
computational description, the user must examine the computer listing for 
CNFWBK in Appendix B, but some general computational details are given in 
Section 6.2. 5.2. A functional flow chart of CNFWBK is shown in Figure 6.4-2. 

The computational flow pattern for CNFWBK is identical to CNBACK 

with the only difference between the routines being the diffusion-node 

temperature finite-difference algorithm. On the first iteration only the 

source locations zeroed out and the present temperatures stored, VARIABLES 1 

is called and variable C. impressed source q. and variable coefficients 

r, ■ X ■ 

6^ (diffusion-diffusion and diffusion-arithmetic) evaluated. All quantities 
which are evaluated at time, t^, ,are summed in accordance with equations 
(6,4-6) and (6.4-8). CSGMIN is evaluated and the diffusion-node tempera- 
tures calculated; note the arithmetic-node temperatures are not calculated 
on the first iteration. 

On the second and succeeding iterations the quantities C^^ q^ and 
(diffusion-diffusion and diffusion-arithmetic) are not updated. 
Diffusion-node temperatures are calculated and DKLXCC determined. Every 
third iteration, if a diffusion-node temperature is converging, a linear 
extrapolation to accelerate convergence is performed (refer to Section 6.2.7). 
If arithmetic nodes are encountered, the appropriate q^ and (for arithme- 
tic nodes) are evaluated once per time-step. Arithmetic-node temperatures 
are calculated and ARLXCC determined. 

Control constants DRLXCC and ARLXCC are checked against DRLXCA and 
ABLXCA, respectively each time-step; if both criteria are satisfied the 
iterations cease, otherwise the iterations continue NL00P times and the 
message "RELAXATION CRITERIA NOT MET" is printed. 

Diffusion-node and arithmetic-node temperature changes between 
time-steps are calculated and stored in DTMPCC and ATMPCC, respectively. 
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If DTMPCC > DTMPCA or if ATMPCC > ATMPCA, the just completed calculations 
are erased and the time-step shortened (refer to Section 6. 2. 5.2). 

6. 4. 2. 4 Control Constants 

The control constants for CNFWBK are used in exactly the same 
way as used in CNBACK. Control constants ARLXCA, CBLXCA, DTIMEI, I^00P, 
0UTPUT, fflid TIMEI5D must be specified as indicated in Table 6.2-5 and as 
- described in Section 6. 2. 3. 2; otherwise "run" will terminate with an 
appropriate error message. Ihe function of optionally specified control 
constants ATMPCA, BACKUP, DAMPA, DAMPD, DTIMEH, DTMPCA, and TIME0 is 
described in Section 6.2. 3.2. 

Specification of time-step DTIMEI should be done in conjtmction 
with control constant HL00P which represents the maximxim manber of com- 
putational iterations during each time-step. Since each iterative calcula- 
tion is essentially equivalent to a time-step calculation for an explicit 
method, the combination of DTIMEI and NL00P for a given time period should 
be less than the total number of time-steps by the explicit method for the 
same time period. Note also that TIME0 may be set negative. Specification 
of ARLXCA and DRLXCA depends upon the problem but a typical value is 0.1. 

6. 4. 2. 5 Error and Other Messages 

If control constants ARLXCA, DRLXCA, DTIMEI, NL00P, 0UTPUT and 
TIMEND are not specified the following error message will be printed for 
each. 


ARLXCA 

"N0 ARLXCA" 

DRLXCA 

"N0 DRLXCA" 

DTIMEI 

"N0 DTIMEI" 

NL00P 

"N0 NL00P" 

0UTPUT 

"N0 0UTPUT INTERVAL" 

TIMEND 

"TRANSIENT TIME N0T SPECIFIED 


If the long pseudo-compute sequence LPCS is not specified, the 
error message will be, 

"CNFWBK REQUIRES L0N6 PSEUD0-C0MPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficient 
(NDIM < (3* NND + NNA + NNB)), the message will be, 

" LOCATIONS AVAILABLE" 
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llote that the number presented will be negative Indicating the additional 
storage locations required. 

If CSGMIN < 0, the following message will be printed, 

"CSGMIN ZER0 OR NEGATIVE" 

If either ARLXCA or DRLXCA is not satisfied with NL00P iterations, 
the following message will be printed, 

"RELAXATI0N CRITERIA N0T MET" 

Checks on the control constants, the pseudo-compute sequence and 
the dynamic storage allocation are made in the following sequence with 
the "run" terminating if a single check is not satisfied, 

NLOOP, TIMEND, OUTPUT, ARLXCA, DTIMEI, DRLXCA, LPCS and 

dynamic storage allocation. 
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Table 6.4-2. 


Basic Computational Steps for CNFIaIBK 


1. Specification of control constants. Control constants ARLXCA (if 
m > 0), DBLXCA (if IIND > 0), DTIMEI, NL00P, 0UTPUT and TIMEND 
(TIMEND ^ TIMEO) must be specified. LPGS is required. (Refer to 
Table 6.2-5 for values and Section 6. 2. 3. 2 for description. 

2. Sufficiency check on d3mainic storage. 

Requirements = 3*NND + NBA + ]SNB (NBD » diffusion nodes, NBA = arithmetic 
nodes and NNB == boundary nodes) 

-3. Setting and/or calculation of time-step, At. (Refer to Section 6.2.4 
for detailed procedure.) Time-step » DTIMEI. 

4. Setting of iterative D0 loop, 1 to NL00P. 

5. Setting of source locations to zero. 

6. Calling of Variables 1. (Refer to Section 6. 2. 2. 2 for description.) 

7. Checking of BACKUP. (Refer to Section 6. 2. 3. 2 for description.) 

8. Diffusion-node temperature calculations, first iteration only. Evalua- 
tion of qi, C^ and Damping of radiation heat transfer. (Refer to 

' Section 6. 2. 5. 2.) Calculation of diffusion-node temperature. The com- 
putational algorithm depends upon the presence of radiation heat trans- 
fer, but the method of solution is the. Crank-Nlcolson algorithm (half 
forward and half backward, refer to Section 6. 2. 5. 2). 

9. Conversion of T^ to “R (Rankine) . 

10. Diffusion-node temperature calculation, successive iterations after 

first. Repeating of step S except that and Gj^ are not updated. 

Calculation of DRLXCC. 

11. Acceleration of convergence every third Iteration if linear extrapola- 
tion is met (refer to Section 6.2.7). 

12. Conversion of T^ to degrees Fahrenheit. 

13. Calculation of arithmetic-node temperatures, second and succeeding 
iterations; arithmetic-node temperatures are not calctilated on the 
first iteration (refer to Section 6, 2, 5.2 for details). 

14. Conversion of temperatures to degrees Rankine. 

15. Checking of ARLXCA and DRLXCA for convergence and 0PEITR for output. 

If both ARLXCA and DRLXCA are satisfied, iterations during a time- 
step cease, otherwise NL00P iterations are performed. 

16. Checking of ATMPCA and DTMPCA. If either one is not satisfied time- 
step is shortened, previous temperatures erased, and temperatures 
recalculated for shortened time-steps (refer to Section 6. 2. 5. 2). 

17. Conversion of temperatures back to degrees Fahrenheit. 

18. Calling of VARIABLES 2 and checking of BACKUP (refer to Section 6. 2. 2. 3 
and 6. 2. 3. 2). 

19. Advancing of time, checking Of time to print, and the printing of the 
output interval. 

20. Calling of 0UTPUT CALLS. 

21. Checking for problem end-time stored in user specified control constant 

TIMEND. ■ 
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6.4.3 Subroutine; CNVAEtB 
6.4. 3.1 General Comments 

Subroutine CNVARB uses an Implicit finite difference algorithm 
that is a composition of forward-differencing and backward-differencing. 
The proportion of forward to backward to be used is calculated internally 
by using a weighting factor, 0, that is dependent upon the ratio of the 
explicit stability criterion as stored in the control constant CSGMIN 
divided by the computational time-step stored in DTIMEU. The weighting 
factor can vary each time-step but is constrained to range, 0 ^ 1/2 

(refer to Section 6. 2. 5. 2 or Section 6. 4. 3. 2). A 3 of one-half yields 
the Crank-Nlcolson half-forward and half-backward expression, whereas a 
3 of zero yields the standard backward-difference expression. 

Except for the weighting factor, 3, the computational procedure 
and the use of the various control constants in CNVAjRB is essentially 
identical to subroutine CNFWBK. 


Solution characteristics should be very similar to CNFWBK 
solutions with expectation that CNVAEB solutions would be more optimum in 
terms of accuracy and solution time. Solutions are not presently available 
to verify or refute the expected advantages of CNVARB solutions. 


6. 4. 3. 2 Finite Difference Approximation and Computational Algorithm 


The numerical solution algorithm used in subroutine CNVARB is 
a combination of forward-differencing and backward-differencing with the 
weighting of each determined by the ratio of control constants CSQIIN/DTIMEU. 


The combination forward-backward differencing with weighting 
can be expressed as: 

^4 ■ / P \ P 

+ tt.o - 6) 

1 * 1,2,. « . ,N 
n “ nth time-step 
3 “ weighting factor (0 < 3 ^ 1/2) 
’^j,n+l “ <^°“stant, N < j < p 
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If equation (6.4-9) is multiplied by 2,0 and the known quantities 
(at time-step, n) and the unknown qtiantities (at time-step, n+1) separated, 
the algorithm used in CNVARB may be obtained readily. 

Diffusion Nodes 


Diffusion-node temperatures are solved by "successive point" 
iteration. The finite difference iterative form as used in CNVARB can 
be found by multiplying equation (6.4-9) by 2-0 and by using appropriate 
time-step, n, and iteration, k subscripts. 


\,k+i ■ 


[Q - (q.) 1/G 

■•^sum ^^x'^ave-' simi 


(6.4-10) 


where, Q. = 2 q. + 2 C. T 


x.n 


i,n i,n 


/ P , P 4 4 \ 

+ eW Z a.. (T. - T, )+ Z Ob.. (T? - TT )) 

xj.n ij,n 3 ,n i,n / 

Q. + (2.0-3’)f Z G.. T. , + Z G.. T. ,\ 

^ J»k+1 


( 6 . 

(6.4- 


* 2 C. + (2.0- g») Z a.. 
Sian i,n xj ,n 


(6.4-13) 


G . . « a . . + Ob . , TT n 

ij ,n ij ,n ij,n 


(6.4-14) 

(£ = k, if j ^ i and f => k+1, if j < i] 


<Vave ■ “‘’y.a "<.k> 


(6.4-15) 


average heat loss from the ith node, called radiation 
damping (refer to Section 6.2.6 for details) 

- 0, if radiation is not present 

B*= 2.0*CSGMIN/DTIMEU (range allowed, 0 <_ g'< 1.0, note 3’ = 2g) 

n = nth time-step; k = kth iteration 

C^,q^,a^j ,b^^ = optionally specified (refer to Tables 6.2-1 - 6.2-4) 

C. * C. /At 
i,n i,n 

i = 1,2,...,NND 

T, ; T. , = constant, (NND + NNA) < j < p (p is the total number of 

j ,n 3 ,K 

nodes and NNA is the number of arithmetic nodes) 
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Arithmetic Nodes 


Arithmetic nodes are calculated in the same manner in all the 
SINDA numerical solution routines. For the finite difference algorithm 
refer to either Section 6. 3.1*2 or Section 6. 2. 5. 2. 

6.4. 3.3 Comments on the Computational Procedure 

The important steps of the computational procedure used in sub- 
routine CNVARB are indicated in Table 6.4-3. For a detailed step-by-step 
computational description, the user must examine the computer listing for 
CNVARB in Appendix B, but some general computational details are given in 
Section 6.2. 5.2. A functional flow chart of CNVARB is shown in Figure 6.4-3. 

The computational flow pattern for CNVARB is very similar to 
CNFWBK or CNBACK; the slight difference is shown in the flow chart of 
Figure 6.4-3. The basic difference between CNVARB and the other two 
implicit routines is the use of a variable beta, g*, which is calculated 
internally by the routine. Thus, the updating of the variable capacitance 
C^, the impressed source and the variable coefficients (a^^ for con- 
duction and for radiation) during the first iteration and the sub- 

sequent calculation of diffusion-node temperatures in subsequent iterations 
are identical to CNFWBK except for the finite difference algorithm. Use 
of the various control constants and checks are identical to CNFWBK. 

6.4. 3.4 Control Constants 

Control constants for CNVARB are used in exactly the same way as 
used in CNFWBK. Control constant ARLXCA, DRLXCA, DTIMEI, NL00P, 0DTPUT, 
and TIMEND must be specified as indicated in Table 6.2-5 and as described 
in Section 6. 2. 3. 2; otherwise ”run" will terminate with an appropriate 
error message. The function of optionally specified control constants 
ATMPCA, BACKUP DAMPA, DAMPD, DTIMEH, DIMPCA and TIME0 is described in 
Section 6.2. 3. 2. 

6. 4. 3. 5 Error and Other Messages 

If control constants ARLXCA, DRLXCA, DTIMEI, NL00P, 0DTPDT and 
TIMEND are not specified, the following error message will be printed for 
each. 
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APiXCA 

"N0 ARLXCA" 

DRLXCA 

"N0 DRLXCA" 

DTIMEI 

"N0 DTIMEI" 

NL00P 

"N0 NL00P" 

0DTPUT 

"N0 OUTPUT INTERVAL" 

TIMEND 

"TRANSIENT TIME N0T SPECIFIED’ 


-If the long pseudo~compute sequence LPCS Is not specified, the 
error message will be, 

"CNVARB REQUIRES L0NG PSEUD0-C0MPUTE SEQUENCE'' 

If the d 3 mamic storage allocatipn is not sufficient 
(NDIM < (3*NND + NNA + NNB)), the error message will be, 

" L0CATIONS AVAILABLE" 

Note that the number presented will be negative indicating the additional 
storage locations required. 

If CSGMIN £ 0, the following message will be printed, 

"CSGMIN ZER0 or NEGATIVE" 

If either ARLXCA or DRLXCA is not satisifed with NL00P iterations, 
the following message will be printed, 

"RELAXATION CRITERIA N0T MET" 

Checks on the control constants, the pseudo-compute sequence 
and the dynamic storage allocation are made in the following sequence 
with the "run" terminating if a single check is not satisfied, 

NL00P, TIMEND,0UTPUT, ARLXCA, LPCS and dynamic storage 
allocation. 
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Table 6.4-3. Basic Comptitational Steps for CNVARB 

1. Specification of control constants. Control constants ARLXCA (if 
NM > 0), DRLXCA (if NND > 0), DTH-IEI, NL00P, 0UTPUT and TIMEND 
(TIMEND ^ TIME0) must be specified. LEGS is required. (Refer to 
Table 6.2-5 for values and Section 6. 2. 3. 2 for description.) 

2 . Sufficiency check on dynamic storage. Requirements = 3*NND + NNA + NNB 
(mO) = diffusion nodes, MA ® arithmetic nodes and NHB * boundary nodes) . 

3. Setting and/or calculation of time-step, At. (Refer to Section 6.2.4 
for detailed procedure.) Time-step * DTIMEI. 

4 . Setting of iterative D0 loop , 1 to NL00P . 

5. Setting of source locations to zero. 

6. Calling of Variables 1 (refer to Section 6. 2. 2. 2 for description). 

7. Checking of BACKUP (refer to Section 6. 2. 3. 2 for description). 

8. Diffusion-node temperature calculations, first iteration only. 

Checking of stable stability criteria. 

Calculation of weighting factor 8'®= 2 , 0*CSGMIN/DTIMEU . (0 £ 3’<. 1*0) 

Conversion of temperatures to degrees Rankine. 

Damping of radiation heat transfer (refer to Section 6. 2. 5. 2). 

Calculation of diffusion-node temperatures using forward-backward 
algorithm with variable beta (8’). 

Calculation of DRLXCC. 

9. Diffusion-node temperature calculations, successive iterations after 
first. Repeating of step 8 except that q^, Ci and Gjj. are not updated. 
Calculation of DRLXCC. 

10. Acceleration of convergence every third iteration if linear extrapolation 
criterion is met (refer to Section 6.2.7). 

11. Conversion of T^ to degrees Fahrenheit, 

12. Calculation of arithmetic-node temperatures every iteration (refer to 
Section 6. 2. 5. 2 for details). 

13. Conversion of temperatures to degrees Rankine. 

14 . Checking of ARLXCA and DRLXCA for convergence and 0PEITER for output . 

If both AEILXCA and DRLXCA are satisfied, Iterations during a time-step 
cease, otherwise NL00P Iterations are performed. 

15. Checking of ATMPCA and DTMPCA. If either one is not satisfied time- 
step is shortened, previous temperatures erased, and temperatures 
recalculated for shortened time-steps (refer to Section 6. 2. 5. 2). 

16. Conversion of temperatures back to degrees Fahrenheit. 

17. Calling of VARIABLES 2 and checking of BACKUP (refer to Section 6. 2. 2. 3 
and 6, 2. 3. 2). 

18. Advancing of time, checking of time to print, and the printing of the 
output interval. 

19. Calling of 0UTPUT CALLS. 

20. Checking for problem end time stored in user specified control constant 
TIMEND. 
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6.5 


Steady State Niimerical Solution Routines 

SXNBA steady state nvunerical solution routines number three. 

These steady state routines are identified as follows: 

CINDSS Block iterative method 

Requires short pseudo-compute sequence (SPCS) 

Clin)SL Successive point iterative method 

Requires long pseudo-compute sequence (LPCS) 

CINDSM Modified CINDSL for radiation -dominated problems 
Requires long pseudo-compute sequence (LPCS) 

A detailed description of steady state routines is presented in 
the pages to follow with liberal reference to materials presented in 
Section 6.2. A brief description of these routines follows. 

CIMDSS which uses the short pseudo-compute sequence (SPCS) was 

the first steady state routine developed for SINDA (via CINDA and CINDA-3G) 

as a result, some of the features contained in subsequent steady state 

routines are not used in CINDSS. If a transient analysis is to be performed 

following a steady state analysis, CINDSS must be used with a transient 

routine that also requires SPCSi The ’’block" iterative method (refer to 

Section 5.2.3) used by CINDSS should lend Itself to some types of problems 

4 4 

which are highly nonlinear with terms such as G. . (T, - T,). With ’’block" 

ij j i 

iteration, both T^ and T^ are changed simultaneously. Solution convergence 
is based upon a temperature relaxation criterion stored in DRLXCA for 
diffusion nodes and ARLXCA for arithmetic nodes. 

CINDSL requires the long pseudo-compute sequence (LPCS) and uses 

the "successive point" iteration method (refer to Section 5.2.3). Any 

transient analysis routine coupled with CINDSL must require LPCS. Solution 

time for CINDSL is less theai CINDSS; as a result, it is used more often 

than CINDSS. A major problem with CINDSL is that a highly nonlinear 

problem can present convergence difficulties unless considerable amount of 

damping is used. For example, a radiation-dominated problem contains many 
4 4 

ab^j (Tj - Tj_)* With "successive point" iteration, T^ may be updated and 
Tj^ not for a given conductor; as a result, the resultant heat flow 
calculation could present difficulties because of large change in values. 
CINDSL has the acceleration of convergence feature, whereas CINDSS does not 
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Solution convergence is based upon temperature relaxation criterion stored 
in DRLXCA for diffusion nodes and AELXCA for arithmetic nodes. 

CINDSM is the latest addition to the SINDA library of steady 
state routines. CINDSM requires the long pseudo-compute sequence and uses 
"successive point" iteration. The routine was specifically developed to 
solve radiation-dominated problems. Solution convergence is based upon 
system energy criterion stored in BALENG. 



6.5.1 


Subroutine: CINDSS 


6. 5. 1.1 General Comments 

Subroutine CINDSS is a steady state routine that requires the 
short pseudo-compute sequence (SPCS) and ignores the capacitance values of 
diffusion nodes to calculate steady state temperatures. Diffusion nodes are 
solved by a '’block*' iterative method as discussed in Section 6. 5. 2. 3, whereas 
arithmetic nodes are solved by a "successive point" iterative method also 
discussed in Section 6. 5. 2. 3. For steady state solutions diffusion nodes are 
not necessai^r; as a matter of fact, solutions will be achieved more quickly 
if all diffusion nodes are specified as arithmetic. The use of diffusion 
nodes in a steady state solution allows for the direct use of the transient 
model. 

A series of steady state solutions at various points in a time 
period can be accomplished by specifying control constants TIMEN and 0UTPUT. 
0UTPUT is used both as the output interval and the computational interval. 

The instiructions with the appropriate call are made in VARIABLES 1 to modify 
boundary conditions with time. 

The CINDSS call can be followed by a call to one of the transient 
solution subroutines which has the same short pseudo-compute sequence 
requirements such as GNFRWD. In this manner the steady state solution 
becomes the initial conditions for the transient analysis. It is important 
to remember that control constants specified for the steady state routine 
will be used by the transient routine unless initialized to the desired 
values. Since CINDSS utilizes control constants TIMEND and 0UTPUT for the 
steady state-transient problem, the user must specify their values in the 
execution block after the steady state call and prior to the transient 
analysis call. CINDSS does not utilize the acceleration of convergence 
feature as discussed in Section 6.2.7. 

Solution convergence is based upon a temperature relaxation 
criterion stored in control constants DRLXCA for diffusion nodes and ARLXCA 
for arithmetic nodes. Normally, identical values are specified for both 
DBLXCA and ARLXCA. Sufficient information is not presently available to 
Indicate different values for DRLXCA and ARLXCA, A method to indicate the 
accuracy of the "converged" temperatures is not presently available. It 
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should also be noted that '’converged" temperatures could have large system 
energy unbalance. 

6. 5. 1.2 Finite Difference Approximation and Computational Algorithm 

The steady state heat balance equation at the ith node may be 
readily expressed as. 


P 

q . + 2 
j-1 




P 

- T.) + 2 

j=l 




i « 1,2, ..,N 

Tj * constant, N < j _< p 


(6.5-1) 


Equation (6.5-1) represents a set of nonlinear algebraic equations 
to be solved simultaneously. Since CINDSS solves temperature of nodes 
specified as diffusion (nodes with capacitance even though a steady state 
solution is desired) by the "block" iteration method and temperatures of 
nodes specified as arithmetic (no capacitance) by the "successive point" 
iteration method, two successive approximation algorithms are used. 


Diffusion Nodes (if any) 


T. T,,i * DD* T. , 
i,k+l l,k 


‘ %.k 

+ l-i . : 


P 

S G 
j=l 


(6.5-2) 




where, k = kth iteration; i = 1,2,..., NND (number of diffusion nodes) 
q^,ay ,b^j = may be optionally specified (refer to Tables 6.2-*l - 6.2-4) 

constant, (NND + NNA) < j <. p (NNA is the nvariber of arithmetic 


j»k 

nodes and p is the total ntimber . of nodes) 

^ 4 - V “ a.. , + Ob.. , (T? ^ + T? ,)(T. . + T. , ) 

i3»k xj,k xj,k ' j,k x,V' j,k x,k' 

DN = DAMPD (diffusion node damping factor) 

DD = 1.0 - DN 


Arithmetic Nodes (if any) 


^1,W1 - “* \,k + 


<%,k + j, “ij.k "3,k+l ^ °«,k "j,k> 


P 

Z G . , , 

3-1 
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where, k = kth iteration; i = (NND -i- l),(mn) + 2),..., (NHD + NNA) 
*^i**ij*^ij ~ optionally specified (refer to Tables 6.2-1 - 6.2-4) 

T , = constant, (NND + NBA.) < J £ p (MA is the number of 

J,K 

arithmetic nodes and p is the total number of nodes 

=ia.k - “ij.k + + •'i.k) 

(Jl » k, if j > i and I = Icfl, if j < i) 
AN = BAMPA (arithmetic node damping factor) 

AD = 1.0 - AN 


6. 5. 1.3 Comments on the Computational Procedure 

The important steps of the computational procedure used in the 
steady state subroutine CINDSS are indicated in Table 6.5-1. For a detailed 
procedural description, the user must examine the computer listing for CINDSS 
in Appendix C, but some general computational details are given in Section 
6. 2. 5. 3. A functional flow chart of CINDSS is shown in Figure 6.5-1. The 
user is required to specify the maximum number of iterations to be per- 
formed via control constant NL00P and the diffusion-node temperature change 
relaxation criteria DRLXCA and the arithmetic-node temperature change 
criteria ARLXCA. The iterations continue \mtil either NL00P is satisfied 
or both DRLXCA and ARLXCA are satisfied. If DRLXCA and ARLXCA are not 
satisfied with NL00P iterations, an appropriate message is printed. 

VARIABLES 1 and 0DTPUT CALLS are performed at the start and VARIABLES 2 and 
0UTPUT CALLS are performed upon completion. Control constants DAMPD for 
diffusion nodes and DAMPA for arithmetic nodes are so-called damping factors 
which are multipliers of the ’’new” temperatures; the factor 1.0 — DAMPD 
(or 1.0 - DAMPA) is a multiplier for the ”old’* temperatures. This 
weighting of "old” and "new" temperatures is useful for damping oscilla- 
tions due to nonlinearities. For nonlinear systems, the damping factors 
are specified to be less than one. If not specified, the damping factor 
is set to 1.0. As a point of interest, it appears that if a linear system 
is to be solved, the convergence could be accelerated by using the damping 
factor greater than one. The diffusion nodes receive a "block" iteration, 
whereas the arithmetic nodes receive a "successive point" iteration; 
acceleration features are not utilized. 
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6.5.1. 4 Control Constants 

)' ■ — 

' Control constant NL00P must be specified and control constants 

ASLXCA and DRLXCA must be specified if NNA > 0 and NND > 0, respectively; 
otherwise "run" will terminate with an appropriate error message. Control 
constants DAMPA and DAMPD may be optionally specified among others. Control 
constant characteristics are tabulated in Table 6.2— 5 and description of 
these control constant is presented in Section 6. 2. 3. 2. Specification of 
NL00P is dependent upon the values of ARLXCA and DRLXCA and thus the 
accuracy of solution. Since the type of problem will influence accuracy, it 
appears that a trial and error procedure is the only practical way of 
determining realistic control constant values. 

6. 5. 1.5 Error and Other Messages 

If control constants ARLXCA, DRLXCA and 15L00P are not specified, 
the following error message will be printed for each, 

ARLXCA "N0 ARLXCA" 

DRLXCA "N0 DRLXCA" 

NL00P "N0 NL00P" 

I ' ' 

If the short pseudo-compute sequence SPCS is not specified, the 
error message will be, 

"CINDSS REQUIRES SH0RT PSEUDO-C0MPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficent (NDIM < NND) 

will be, 

" L0CATIONS AVAILABLE" 

Note that the number printed will be negative indicating the additional 
storage locations required. 

If both temperature change relaxation criteria ARLXCA and DRLXCA 
are not met with NL00P iterations, the message will be, 

"ITERATI0N C0UNT EXCEEDED, L00PCT = _" 

Checks on the control constants, the pseudo-compute sequence, 
and the dynamic storage allocation are made in the following order with the 
"run" terminating if a single check is not satisfied. 

NL00P, ARLXCA, DRLXCA, SPCS, and dynamic storage allocation. 
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1 . 

2. 

3. 


4. 

5. 

6 . 
7. 


8 . 

9. 


10 . 

11 . 


12 . 

13. 

14. 


Table 6.5.1. Basic Computational Steps for CIHDSS 

Specification of control ccaistants. Control constants AMiXCA (if KM > 0), 
BRIaXCA (if NND > 0) and NL00P must be specified. SPCS is required. 

(Refer to Table 6.2-5 for values and Section 6. 2. 3. 2 for description.) 

Sufficiency check on dynamic storage. Requirements = NHD (KND = 
diffusion nodes). 

Setting of TIMEN for first iteration and succeeding iterations. 

TIMEN = TIME0, first iteration 

TIHEN = TIME0 + 0UTPIJT, succeeding iterations 

Setting of iterative loop for all nodes, kl = 1, NL00P 

Setting of source locations to zero. 

Calling of VARIABLES 1 (refer to Section 6. 2. 2.2 for description). 

Calculation of diffusion-node temperatures by "block” iteration if 
NND > 0 (refer to sections 6. 2. 5. 3 and 6. 5.1. 2). 


’^l.k+1 - \.k ^ 


P 

Z G 

3»1 


ij,k 


DN = DAMPD and DD * 1.0 - DN 


Calculation of DRLXCC. 


Calculation of arithmetic-node temperatures by "successive point" 
iteration if NM > 0 (refer to Sections 6. 2. 5. 3 and 6. 5. 1.2). 


T * AD* T + 

^l,k+l ^ ^i,k 


AN* (q^ ^ + Z G,.^ ^ T 


i.k ’4:, "j»k+l ■ 4-44.1 ''ij»k ^j,k 


+ Z G,, ^ T, J 


P 

Z G 

j=l 


i3,k 


AN = DAMPA 

AD = 1.0 - DAMPA 


Calculation of ARLXCC. 

Checking of DRLXCC and ARLXCC against the relaxation criteria DRLXGA 
and ARLXCA, respectively, for convergence. If both ARLXGA and DRLXCA 
are satisfied, iterations cease, otherwise NL00P iterations are 
performed. 

Calculation of system energy balance which is stored in ENGBAL. 

Call VARIABLES 2 and 0DTCAL, print ENGBAL and L00PCT. 

Check if TIMEND = TIMEN. 
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CIHDSS 


TIMEN + TIME0 (1st time) 

TIMEN = TIME0 + 0UTPUT (succeeding) 


Iterative X00P 
K1 = 1,NL00P 


Set all sources to zero 


Call VARIABLES 1 


Call OUTPUT CALLS, first pass only 


Update Gj^ (diffusion-diffusion, diffusion- 
arithmetic) ; compute diffusion-node temperatures 
by block iteration method 


Update G-^ (arith.-arith. )» compute arith.- 
node temperatures by successive point iteration 


Converge 


Call 0UTPUT CALLS 
if 0PEITR 0 0.0 


K1 = HL00P 


Print message 


Calculate energy balance 


Call VARIABLES 2 


TIME0 = TDIEN 


Call 0UTPUT CALLS 


Return] 


Figure 6.5-1. Functional Flow Chart for CINDSS 


















6.5.2 


Subroutine; GINDSL 


6. 5. 2.1 General Comments 

Subroutine CINDSL is a steady state routine that requires the 
long pseudo-Gompute sequence (LPCS), Both diffusion- and arithmetic-node 
temperatures are calculated by a "successive point” iteration computational 
technique. Every third iteration a linear extrapolation is performed to 
accelerate convergence. CINDSL generally yields significantly faster 
solutions than CINDSS, but nonlinear problems such as those with radiation 
heat transfer can pose considerable convergence difficulties unless a large 
amount of damping (low values of DAMPA and DAMPD) is imposed. 

A series of steady state solutions at various points in time can 
be generated by specifying control constants TIMEND and 0UTPUT. 0UTPUT 
is used both as the output interval and the computation interval; this 
requires appropriate ralla in VARIABLES 1 to modify boundary conditions 
with time. 

CINDSL can be followed by a call to one of the transient numerical 
solution routines which have the same LPCS requirements. Used in this 
manner the steady state solutions become the initial conditions for the 
transient analysis. Note that since CINDSL utilizes control constants 
TIMEND and 0UTPUT for the coupled steady state-transient problem, the user 
must specify the values of TIMEND and 0UTPUT in the execution block after 
the steady state call and prior to the transient analysis call . 

Solution convergence is based upon a temperature relaxation 
criterion stored in control constants DRLXCA for diffusion nodes and ARLXCA 
for arithmetic nodes. Normally, identical values are specified for both 
DRLXCA and ARLXCA for lack of anything better. The damping factors DAMPD 
for diffusion nodes and DAMPA for arithmetic nodes are merely multipliers 
of "new" temperatures and the factor 1.0 - DAMPD (or 1.0 - DAMPA) is a 
multiplier of the "old" temperatures. Normally, these damping factors are 
specified to be less than 1.0, but for a linear system the convergence 
probably could be accelerated by using a damping factor greater than one. 

6. 5. 2. 2 Finite Difference Approximation and Computational Algorithm 

The set of steady state heat balance equations, 
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i « 1,2,..., K 


Tj ~ constant N < j <^ p 


is solved by a re-iterative scheme called a "successive point*' iterative 
method here. Both diffusion-node and arithmetic-node temperatures are 
solved in this manner. The only difference between the two algorithms is 
that control constant DAMPD is used with diffusion nodes and control 
constat DAMPA is used with arithmetic nodes. 

Diffusion Nodes (if any) 




A ,= 1+1 


P 

E G 


ij,k 


where, 1 = 1,2,...,NHD; k = kth iteration j 

q^,a^^,b^^ * may be optionally specified (refer to Tables 6.2-1 - 6.2-4) 

T. . = constant, (NMD + NNA) < j < p (NNA is the number of arithmetic 
J f ^ 

nodes and p is the total number of nodes) 

DN * DAMPD (diffusion-node damping factor) 

DD = 1.0 - DN 

e«,k " + "i.k> 


(£ = k if j > i and Jl = k+1 if j < i) 


Arithmetic Nodes (if any) 


\,kfl = ^i,k + 


^^i,k ®ij,k ^j,k+l ‘^ij,k '^j,k^ 


P 

E 6 


ij,k 


where, i = (NND + 1), (NND + 2), . . ., (NND + NNA) 

*^i**ij*^ij* optionally specified (refer to Tables 6.2-1 - 6.2-4) 

k * constant (NND + NNA) < j ^ p (NNA is the number of arithmetic 
nodes and p is the total number of nodes) 

AN “ DAMPA (arithmetic-node damping factor) 

AD = 1.0 - AN 

®ij,k “ ^ij.k ®^ij,k ^ \,k^ 

(£ = k if j > i and Jl * k+1 if j < i) 
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5.2.3 Comments on the Computational Procedure 

The Important steps of the computational procedure used in 
the steady state subroutine CIKBSL are indicated in Table 6.5-2. For a 
detailed procedural description, the user must examine the computer 
listing for CINDSL in Appendix C, but some general computational details 
are given in Section 6. 2. 5. 3. A functional flow chart of CI]:^SL is shown 
in Figure 6.5-2. 

The computational pattern of CINDSL is very similar to CINDSS 
with the differences being that CINDSL uses the long pseudo-compute 
sequence, whereas CINDSS uses the short pseudo-compute sequence, and that 
CINDSL contains the acceleration convergence feature, whereas CINDSS does 
not. The user is required to specify the maximum ntanber of iterations 
to be performed via control constant NL00P and the diffusion-node tempera- 
ture change relaxation criteria DRLXCA and the arithmetic-node temperature 
change relaxation criteria ARL3CCA. The iterations continue until either 
NL00P is satisfied or both DRLXCA and ARLXCA are satisfied. If DRLXCA and 
ARLXCA are not satisfied with NL00P, an appropriate message is printed. 
Acceleration of convergence is performed every third iteration if a 
temperature is converging over two time-steps. 

6.5. 2. 4 control Constants 

Control constant NL00P must be specified and control constants 
ARLXCA and DRLXCA must be specified if NNA > 0 and NND >0, respectively; 
otherwise ’'run" will terminate with an appropriate error message. Control 
constants DAMPA and DAMPD may be optionally specified among others. Con- 
trol constant characteristics are tabulated in Table 6.2-5 and description 
of these control constants is presented in Section 6. 2. 3. 2. Specification 

of NL00P is dependent upon the values of ARLXCA and DRLXCA and thus the 

\ 

accuracy of the solution. Since the type of problem will influence 
accuracy, it appears that a trial and error procedure is the only practical 
way of determining realistic control constant values. 

6. 5.2.5 Error and Other Messages 

If control constants ARLXCA, DRLXCA and NL00P are not specified, 
the following error message will be printed for each. 
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AKLXCA 

DRLXCA 

NL00P 



•’N0 AELXCA" 
"N0 BELXCA” 
"N0 NL00P" 


If the long pseudo-compute sequence LPCS is not specified, 
the error message will be, 

"CINDSL REQUIRES L0NG PSEUD0-COMPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficient, 

(NDIM < 2* (NNA + NND) ) , the message will be, 

" L0CATIONS AVAILABLE" 


Note that the number printed will be negative indicating the additional 
storage locations required. 

"L00PCT = and ENGBAL - " 

If both temperature change relaxation criteria, ARLXCA and 
DRLXCA, are. not met with NL00P iterations, the message will be. 


) 


"ITERATI0N C0UNT EXCEEDED, L00PCT = 

Checks on the control constants, the pseudo-compute sequence, 
and the dynamic storage allocation are made in the following order with 
the "run" terminating if a single check is not satisfied. 


NL00P, ARLXCA, DRLXCA, LPCS, and dynamic Storage allocation. 


\ 
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1 . 


2 . 


3. 


4. 

5. 

6 . 

1 , 


8 . 

9. 


10 . 

11 . 


12 . 

13. 

14. 

15. 


Table 6.5.2 Basic Computational Steps for CIKDSL 


Specification of control constants. Control constants ABLXCA (if 
NNA > 0), DRLXCA (if NND > 0) and NL00P must be specified. LPCS is 
required. (Refer to Table 6.2~5 for values and Section 6. 2. 3. 2 for 
dexcription. ) 

Sufficiency check on dynamic storage. Requirements *= 2* (NND + NNA) 
(NND * diffusion nodes and NNA = arithmetic nodes). 


Setting of TIMEN for first and succeeding iterations. 

TIMEN = T1ME0, first iteration 

TIMEN * TIME0 + 0UTRBT, succeeding iterations 

Setting of iterative loop for all nodes, kl = 1, NL00P. 

Setting of source locations to zero. 

Calling of VARIABLES 1 (refer to Section 6. 2. 2. 2 for description). 

Calculation of dif fusion-node temperatures by '•block" iteration if 
NND > 0 (refer to Section 6.2. 5. 2 and 6.5. 1.2). 


‘i,k+l 


DD* T, , + 
x,k 




p 

Z G.. , T. , ) 
=i+l 


P 

Z 


’ij,k 


DN » DAMPD and DD = 1.0 - DN 


Calculation of DRLXCC. 

Calculation of arithmetic-node teinperatures by "successive point" 
iteration if NNA > 0 (refer to Sections 6. 2. 5. 3 and 6. 5. 1.2). 


T. , * AD* T. , + 

i,k+-l i,k 


"•* («i.k * \i.k =«.k L.k> 




P 

Z G 


ij.k 


Calculation of ARLXCC. 

Checking of DRLXCC and ARLXCC against the relaxation criteria DRLXCA 
and ARLXCA, respectively, for convergence. If both ARLXCA and DRLXCA 
are satisfied, iterations cease, otherwise NLOOP iterations are 
performed. 

Acceleration of convergence each third iteration, if linear extrapola- 
tion criterion is met (refer to Section 6.2.7). 

Calculation of system energy balance which is stored in ENGBAL. 

Call VARIABLES 2 and 0UTCAL, print ENGBAL and L00PCT. 

Check if TIMEND « TIMEN. 


6 - 111 




TIMES = TIME0 + 0UTPUT (succeeding) 


Iterative L00P 
K1 = 1, KL00P 


Set all sources to zero 


Call VARIABLES 1 


Call 0UTPUT CALLS, first pass only 


Update Gjj. (diffusion-diffusion, diffusion- 
arithmetic) j compute diffusion-node temperatures 
by successive point iteration 


Update q^, (arith.-arith.): compute arith.- 
node temperatures by successive point iteration 


Convergec 


Call 0UTPUT CALLS 
0PEITR # 0.0 


= NL00P' 


Print message 


Calculate energy balance 


Call VARIABLES 2 


TIME0 = TIMEN 


Call 0UTPUT CALLS 


T1MENDV ».^ No 
* TIMEN./“ 


Figure 6.5-2 Functional Flow Chart for CINDSL 
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6.5.3 


Subroutine; CINDSM 
6. 5. 3.1 General Comments 

Subroutine CINDSM is a steady state routine specifically generated 
for radiation dominated problems. CINDSM requires the long pseudo-compute 
sequence (LPCS) and is considerably different from CINDSL. CIM)SM is 
based on the use of pseudo linear equations which are the result of 
linearizing the radiation conductors. These equations are solved by using 
the "successive point" method with LAXFAC iterations. Updating of the 
properties as well as the linearized conductors occur outside of the 
iterative loops. Temperature convergence is based on a criterion that is 
continually tight aied until either the NL00P iterations or the system 
energy balance criterion stored in BALENG has been satisfied. 

The acceleration of convergence by linear extrapolation as used 
in CINDSM is essentially the same as used in the other SINDA numerical 
solution routines, but in lieu of limiting the extrapolation by an allowable 
slope value (refer to Section 6.2.7) the maximum tanperature change of the 
network on the last iteration is used as the allowable value. 

Information available at this time Indicates that each problem 
appears to have an optimiua combination of NL00P, DAMPD, and LAXFAC values. 

An NL00P of 100, a DAMPD of 0.5 and a LAXFAC of 10 has been successfully 
applied to spacecraft problems with radiation domination, but the solution 
time is rather long. 

6.5. 3.2 Finite Difference Approximation and Computational Algorithm 

The set of steady state heat balance equations. 


q. + Z a.. (T. -T^ + Z ab , (T^ - T^) = 0 
i i3 3 1 13 3 

i=l,2,...,N 
T^ “ constant, N < j ;< p 

is solved by a re-iterative "successive point" method after linearization. 

4 4 

Linearization is achieved by letting ~ with 

6^ « Cb^j (T^ + T^KTj + T^. This yields 
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(6.5-6) 


P P 

q, + Z a,, (T, - T.) + Z (T, - T.) « 0 
X 13 3 i r 3 i 

Diffusion and Arithmetic Nodes 


No distinction is made between diffusion and arithmetic nodes. 
As a result, the following algorithm applies to both types of nodes. 


\,k “ ■'i.k 


'^J.k+1 * ^J.k' 

'P " 

2 ®ii L 

j=l 


(6.5-7) 


where, i = 1,2, . . . , (NND + NNA); p = total number of nodes 
k = kth iteration 

L = before each LAXFAC iterative loop 
Tj * constant, (NND + NNA) < j i P 

DN * DAMPD (diffusion-nOde damping factor; DAMPA is not used) 

DD = 1.0 - DA]MPD 

>=U.L " + "i,L> 

(G.. is updated once before each LAXFAC iterative loop 
13 »L 

NNA * number of arithmetic nodes 
NND == number of diffusion nodes 

q^,ay,by = may be optionally specified (refer to Tables 6.2-1 - 6.2-4) 
6. 5.3. 3 Comments on the Computational Procedure 


A detailed step-by-step computational procedure as used iii the 
steady state routine CINDSM is presented in Table 6.5-3. For a more 
detailed procedural description, the user must examine the computer 
listing in Appendix C. A functional flow chart that is compatible with 
the step-by-step description of Table 6.5-3 is shown in Figure 6.5-3. 

CINDSM is considerably different from either CINDSS or CINDSL 
because of the use of a variable convergence criterion which is internally 
updated. Overall, from a total system basis, control constants NL00P and 
BALEN6 are the ultimate criteria. 


It should be particularly noted here that unlike CINDSS or 
CINDSL, which use both DAMPA and DAMPD, CINDSM uses only DAMPD. The 
reason for this is that CINDSM does not treat the nodal types as diffusion 
or arithmetic. 
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6. 5* 3. 4 Control Constant 


Control constants BALING, LAXFAC AND NL00P must be specified; 
otherwise the "run" will terminate with an appropriate error message* 

Control constant DAMPD may be optionally specified among others. Control 
constant characteristics are tabulated in Table 6.2,5 and description of 
these control constants is presented in Section 6. 2. 3. 2. Specification of 
BALENG, LAXFAC and NLOOP appears to be a trial and error procedure. 

6. 5. 3. 5 Errot and Other Messages 

If control constants BALENG, LAXFAC, and NL00P are not specified, 
the following error message will be printed for each^ 

BALENG "N0 BALENG" 

LAXFAC "N0 LAXFAC" 

NL00P "N0 NL00P" 

If the long pseudo-compute sequence LPCS is not specified, 
the error message will be, 

"CINDSM REQUIRES L0NG PSEUD0-C0MPUTE SEQUENCE" 

If the dynamic storage allocation is not sufficient, 

(NDIM < ( 3 * NNA + 3* NHD + NGT)), the message will be, 

" LOCATIONS AVAILABLE" 

Note. that the number printed will be nega tive indicating the additional 
storage locations required. 

If either NL00P iterations has been made or if ENGBAL < BALENG, the 
following message is printed, 

"L00PCT = j and ENGBAL " 

Checks on the control constants, the pseudo-compute sequence 
and the dynamic storage allocation are made in the following order with 
the "run" terminating if a single check is not satisfied, 

NL00P, LPCS, BALENG, LAXFAC and dynamic storage allocation. 
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Table 6.5.3. 


Basic Computational Steps for CIKDSM 

1. Specification of control constants. Control constants BALENG, LAXFAC 
and NL00P must be specified. The long pseudo-compute sequence (LPCS) 
is required. (Refer to Table 6.2-5 for values and Section 6. 2.3. 2 for 
description.) 

2. Sufficiency check on dynamic storage. Requirement = 3* (NND + NNA) + NGT 
(NND = diffusion nodes, NNA = arithmetic nodes and NGT = total number 

of conductors). 

3. Setting of TIMEN for the first and succeeding iterations. 

TIMEN = TIME0, first iteration 

TIMEH = TIME0 + 0UTPUT, succeeding . iterations 


4. Constants used in CINDSM 


5. 


6 . 


NLAX = NL00P/LAXFAC (both NL00P and LAXFAC are specified by 
the user) 

RELAX = .05 (initial value used in CINDSM as the allowable 
temperature change) 

DELXXX = .05/NLAX (a nximber used in reducing RELAX for a tighter 
criterion) 

XXXDDM = .001 (a value of RELAX used in CINDSM for a tighter 
criterion) 

* .001/5 (a subsequent value of RELAX for a tighter 
criterion) 

DAMP *= DAMPD (damping factor for all nodes; DAMPA is not used) 


Updating of variables and linearization of radiation . 


Variable q^ and G^ are evaluated by calling subroutine N0NLIN. 

Linearization means that the radiation exchange expressed as crb^. 
Normally, Gij would be updated each iteration as done in CINDSS ^ 
or CINDSL, but in CINDSM G^j is not updated within the D0-L00P 
(kl = 1, LAXFAC) but is updated outside of the loop. 


(Tj-4). 


Iterative D0-L00P (kl = 1, LAXFAC) is established. 

Temperatures of all nodes are calculated by "successive point" iteration 
with no damping . 


A ^j,k 

^i,k+l " p 

j»l 


(6.5-5) 


where, G^j = ai^ + Obij (t| + t|) (Tj + T^) (q^ and Gj^j are not updated 
during the LAXFAC iterations) 

Check on temperature convergence. Temperatures have converged if, 

1^4 ,,.1 - T, , I < RELAX (= .05) 

' i,k+l i,k'max — 


If temperatures have converged, the computation goes out of the 
iteration loop to step (7). 
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Table 6.5.3. (continued) 


Every third iteration, acceleration of convergence is attempted if 
linear extrapolation criterion is met (refer to Section 6.2.7). 

Iteration ceases if LAXFAC Iterations have been performed or if the 
temperatures have converged. 


7. Check on NLAX iterations. 

If in step (6) the number of iterations, LOOPCT > HLAX, the computa- 
tional procedures go to step (9). However, in step (6) if the number 
of iterations LOOPCT < NLAX, then a set of temperature calculations 
is made using "successive point" method with a damping factor and no 
iterations. ^ 

^I,k« ■ \,k ^ ^ i 

Z G.. 
j=l 

where, DN * DAMPD (diffusion node damping factor; note DAMPA is not used) 
= constant 

Allowable temperature change criterion RELAX is reduced to. 


RELAX = .05 - (.05/NLAX) 
and computational procedure goes to step (5). 

8. Repetition of steps (5) through (7) except for temperature convergence 
criterion. 

Temperatures have converged if, 

fcfl " ^i k^ - “ .05/NLAX) 

* max 

9. Assuming step (7) has been satisfied, L00PCT is checked against NL00P. 
If LOOPCT ^ NLOOP, the computation proceeds to step (12) . 

If LOOPCT < NLOOP computation proceeds to step (10). 

10. Reduce RELAX to .001. 

11. Check on temperature convergence. 

If ]T^ - T^ j^l < RELAX (=*= .001) go to step (12). 

jTi fcl ^ X- .001) , LAXFAC is reduced to 

LAXFAC = NL00P - L00PCT, 
and steps p) through (11) are repeated. 

12. Compute system energy balance and store in control constant ENGBAL. 

13. If L00PCT > LAXFAC (original user input value) , go to step (15) 

14. If LOOPCT LAXFAC (original user input value), ENGBAL is checked 
against BALENG. 

If ENGBAL £ BALENG, go to step (14) 

If ENGBAL > BALENG, RELAX is set to, RELAX « .001/5, and steps (5) 
through (14) are repeated with the new RELAX values. 

15. Print ENGBAL; call VARIABLES 2; call 0UTCAL; check if TIMEND » TIME0. 
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6.5-3. Functional Flow Chart for CINDSM 

/ 
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00460 85 * 3000 CONTINUE 

0046 X 85 * END 
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